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Abstract 



The pure gauge theory in 2+1 dimensions is explored, through both a phenomenological 
model and a lattice calculation. The Isgur-Paton model is extended to include a curvature 
term and various mixing mechanisms. The method of inferential statistics is used to 
extract the parameters of best fit and to compare the likelihoods of the various models 
when compared to existing lattice data. The conventional assignment of spin to the 
pscudoscalar state is called into question by the proximity of a spin 4 state in the model, 
which motivates calculating the mass of the spin 4 state on the lattice. Novel lattice 
operators are constructed from a matrix of effective Greens functions which attempt to 
overcome the lattice rotational ambiguities. Correlation functions arc presented for the 
channels with even J, and effective masses extracted. The resulting masses compare well 
with the extended Isgur-Paton model. 
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Chapter 1 



Introduction 



Among the many phenomena currently being explored in the world's particle accelerators, 
few are as interesting as the search for the glueball. Several candidates have presented 
themselves [2] at around 1.6 GeV, but mixing with nearby quarkonia states [3] or other 
non-perturbative effects has kept the situation as clear as mud. In order to identify the 
glueball from the potential candidates, a clear understanding of QCD's nonperturbative 
spectrum is required. Unfortunately, solving QCD completely is still a long ways off. Much 
progress has been made, however, using various techniques. Among the most promising is 
the study of nonperturbative QCD using lattice regularization [4] . Once theory has made 
an unequivocal prediction for the QCD spectrum, experimentalists will be able to verify 
the existence of the glueball. 

While some progress has been made in studying dynamical quarks on the lattice [5] , the 
pure gauge theory without quarks is much better under control and sufhciently interesting 
to warrant study. Without quarks around, the only bound states that exist must be purely 
gluonic: glueballs [6], gluelumps [7,8], and perhaps even more exotic states [9]. Lattice 
methods have produced remarkably accurate estimates of the glueball spectrum [10,11], as 
well as estimates of the color-electric flux proflle [12,13]. Unfortunately, the best estimate 
of the lightest scalar glueball's mass [10], about 1.6 GeV, is smack between the two most 
promising experimental candidates [2]. Thus a more complete understanding of gluonic 
dynamics needs to be developed, so that we may build a model for glueballs which makes 
experimentally testable predictions. 

Several models for glueballs currently exist, each with its own strengths and weaknesses. 
The bag model [14] receives much attention in the literature [7,15,16], but unfortunately 
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is ignorant of the consequences of the no-hair theorem, namely that a spherical shell of 
nonvanishing flux cannot exist [17]. Constituent gluon models also exist [18], but their 
development is hampered by the intrinsically nonperturbative aspect of QCD. The most 
promising candidate for a model of glueballs is the flux tube model, first applied to QCD 
by Isgur and Paton [6] . Chapter [2] reviews the traditional Isgur-Paton model as it applies 
to D=2+l QCD. In this thesis, we work with one less spatial dimension for simplicity - 
there is still plenty of structure to be quite interesting. As is often the case, a system 
reduced to 2+1 dimensions behaves as though it were the cross section of the full 3+1 
system [1], and so understanding glueballs in 2+1 dimensions should impart much insight 
into the full theory. 

Proceeding to Chapter[ni we modify the traditional Isgur-Paton Hamiltonian to include 
a contribution from a curvature term. This term would arise from a flnite coefficient of elas- 
ticity for the flux tube, as is the case with flux tubes in a type-II superconductor [19-22]. 
The traditional Isgur-Paton model does not include states with negative charge conjuga- 
tion eigenvalue in 2+1 dimensions, and so in Chapter [5] we explore several mechanisms 
which give rise to mixing between oppositely orientated flux tubes to produce both the 
positive and negative charge conjugate sectors of the spectrum. Chapter [5] presents the 
results of the calculations, and Chapter [S] explores the behavior of the parameters and the 
merit function for the various models. 

One feature of the flux tube model is the presence of a divergent term in the potential 
for small radius [23]. Taking inspiration from a recent study of flux rings in the dual 
Ginzburg-Landau model [24], Chapter [7] explores how including an effective string tension 
which goes to zero at small radius modifies the spectrum. 

When comparing the Isgur-Paton model's spectrum to existing lattice data [10], some 
surprising features are observed. The model seems to imply that the conventional assign- 
ment of spin to certain lattice operators may be in error. Traditional lattice operators 
belong to one of the symmetry channels of the cubic lattice, which may be compatible 
with several values of spin modulo 4. Usually the lowest value of spin is assigned to the 
lightest state in the channel, but the presence of spin 4 states in the extended Isgur-Paton 
model at the mass of the pseudoscalar {J^'~'' = ^) from the lattice casts doubt on this 
particular assignment, and so in Chapter [5] we begin developing an alternative construc- 
tion of lattice operators that should be free of spin ambiguities. Chapter [TO] details the 
calculation and presents a new algorithm for inverting an Hermitian matrix with the spe- 
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cific band structure encountered here. The resuhs using this new method are presented 
in Chapter [TT] 

We conclude this thesis with Chapter [T^ reviewing the modifications to the traditional 
Isgur-Paton model and their implications for the glueball spectrum and comparing the 
results with the new lattice operators. In Appendix El we present an overview of inferential 
statistics as used in this work. Appendix |B] reviews the details of the Isgur-Paton model 
and demonstrates how states are identified within its context. Finally, representative 
copies of the programs used may be found in Appendix [Cl 

For future reference, here we present the lattice data from [10] in Table [TTTl 





state 


SU(2) 


SU(3) 


SU(4) 


SU(5) 


SU(C50) 


0++ 


4.718(43) 


4.329(41) 


4.236(50) 


4.184(55) 


4.065(55) 


0++* 


6.83(10) 


6.52(9) 


6.38(13) 


6.20(13) 


6.18(13) 


Q+ + ** 


8.15(15) 


8.23(17) 


8.05(22) 


7.85(22) 


7.99(22) 


0-+ 


9.95(32) 


9.30(25) 


9.31(28) 


9.19(29) 


9.02(30) 


2++ 


7.82(14) 


7.13(12) 


7.15(13) 


7.19(20) 


6.88(16) 


2++* 






8.51(20) 


8.59(18) 




2-+ 


7.86(14) 


7.36(11) 


6.86(18) 


7.18(16) 


6.89(21) 


2-+* 




8.80(20) 


8.75(28) 


8.67(24) 


8.62(38) 


1++ 


10.42(34) 


10.22(24) 


9.91(36) 


10.26(50) 


9.98(25) 


1-+ 


11.13(42) 


10.19(27) 


10.85(55) 


10.28(34) 


10.06(40) 


0— 




6.48(9) 


6.271(95) 


6.03(18) 


5.91(25) 






8.15(16) 


7.86(20) 


7.87(25) 


7.63(37) 


Q— ** 




9.81(26) 


9.21(30) 


9.51(41) 


8.96(65) 


0+- 




10.52(28) 


10.35(50) 


9.43(75) 


9.47(116) 


2— 




8.75(17) 


8.22(32) 


8.24(21) 


7.89(35) 


2 — * 




10.31(27) 


9.91(41) 


9.79(45) 


9.46(66) 


2+- 




8.38(21) 


8.33(25) 


8.02(40) 


8.04(50) 


2+-* 




10.51(30) 


10.64(60) 


9.97(55) 


9.97(91) 






9.86(23) 


9.50(35) 


9.65(40) 


9.36(60) 


1+- 




10.41(36) 


9.70(45) 


9.93(44) 


9.43(75) 



Table 1.1: Glueball masses from the lattice in units of the string tension. 



Chapter 2 

The Isgur-Paton Model of 
Glueballs 

The traditional Isgur-Paton model stems from the strong coupling limit of (lattice) QCD [6]. 
In this scenario, the color-electric flux between two point sources does not spread uniformly 
over space but is restricted to the region of a flux tube between the sources. As remarked 
earlier, such a flux tube picture also comes about from the dual Meissncr effect in a 
Ginzburg-Landau type superconductor [19,20,22]. Much numerical evidence [8,13,25] 
points to the existence of this flux tube, even for moderate values of the coupling. The 
flux tube is quantized between the static quarks upto some small cutoff using an adiabatic 
Born-Oppenheimer type approximation. In the limit of infinite quark mass, the model 
provides a good reckoning for mesons [6,8]. 

Gauge invariant bound states also exist in the gauge theory without constituent quarks. 
The simplest topologically is a single loop of color fiux. If the q and q of the flux tube meson 
are brought sufficiently close as to annihilate, the remaining loop of flux might remain as 
a purely gluonic bound state. This closed flux tube is going to have some finite width 
on the order of -^/a, where a is the string tension in units of energy/length. Numerical 
estimates of this width [12] and the distribution of the chromoelectric flux within [25] exist 
and agree with the rough estimate. To simplify matters, though, we will suppose the flux 
tube to be described by a loop of string of vanishing thickness with coordinates given in 
polar form as 

r'=r'(r,$) (2.1) 



A 
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where $ represents the angular variable(s) for either 2 or 3 spatial dimensions. This loop 



is going to have complicated dynamics akin to the Nambu-Goto string [26]: 



'S'Nambu— Goto (X a d Ci/g 



(2.2) 



But we desire a simpler model with more intuitive dynamics. Without the static quarks 
around to provide a natural adiabatic limit, we are a little more cavalier in performing 
the Born-Oppenheimer approximation. Nevertheless, we make the approximation, hoping, 
despite the evidence of previous calculations [27], that the magnitude of the resulting con- 
tributions will justify the gross approximation that is made. The fluctuations of the string 
are divided into an overall radial "breathing" mode of low conjugate momentum and trans- 
verse vibrations of higher momentum. In applying the semi-classical approach to quantiz- 
ing this system, we separate off the presumably "fast" transverse vibrations and quantize 
these as phonons on the string providing a radial potential Vphonon i'r) = Vphonon (f) within 
which the loop fluctuates in its breathing mode. (Details of this procedure can be found 
in [28].) Also arising from the angular quantization of the loop is the Liischer universal 
correction term 



which depends on the number of spatial dimensions D [29]. Evidence for the actual 
existence of this contribution (in the case of a straight string within a meson) has been 
presented in [30]. When we write down the expressions for the energy on our way to 
quantization, we need to include the self-energy of the string. 



which is easily seen to be the product of the string tension and the total (average) length 
of the string. 

This string of flux is related to the flux-ring solution of the dual Ginzburg-Landau 
(dGL) theory [22]. In that theory, the non-perturbative QCD vacuum is represented by 
a type-II superconductor, and vortex field solutions form around a line on which the 
monopole field vanishes. Having no kinetic term, the dGL flux-ring collapses to its center 
and thus cannot represent a stable bound state. Introducing a kinetic term for the flux 

^The reason for the nought on 70 is that later 7 will represent a free parameter with a constant 
correction of 70 




(2.3) 



r 



fj, — 2'Kar 



(2.4) 
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ring meets with varying degrees of success at solving the relativistic Schroedinger equation 
that results. Progress has been made recently, and soon perhaps the detailed dynamics of 
the dGL flux-ring can be matched with the gross features of the flux-tube model [24] . 

The IP model approaches the kinetic term in a non-relativistic way. Consequently, we 
can expect the model's predictions to be off by some amount due to "relativistic effects." 
This systematic error will be later handled by increasing the errors on the "experimental" 
data points taken from the lattice by 5% in quadrature with the statistical errors, so that 
we may say that the model "fits within 10%" of the lattice data. Details will be provided 
later of the effect of this fudge factor. 

For a classical bosonic string [31] of circular radius r and mass fi, we write its kinetic 
energy as 

KE = P^/2iJ, (2.5) 
where Pr is the momentum in the f direction. Substituting the expression for /x, we have 

KE = P^/A-nar. (2.6) 



This operator, as written, will not be Hermitian upon quantization (because it is not 
symmetric). Following standard techniques [32], the independent variable becomes 



(2.7) 



so that, when P^. ihS- and h is set equal to unity. 



-1 P2 

KE = (2.8) 



Attu r 
477(7 

T 3 1/2^ 3 1/2^ 
-t 2^^ dp 21^ dp 

477(7 pV3 pl/3 

-1 9 

4^49^2' 



(2.9) 
(2.10) 
(2.11) 



which is Hermitian. When normalizing the wavefunction, we must take care to include 
the appropriate measure induced by the change of variables: 



5gf = I dp~g{p)f{p) = ^-j r'/^drg{r)f{r). 



(2.12) 
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The quantization of the transverse vibrations into phonons is easily visuaUzed. Vibra- 
tions are classified according to their number of nodes m and their polarization, left-going 
or right-going. We write the occupation number of a particular mode as G Z+. Each 
mode contributes a term ~ ^(n^ + n^)/r to the potential, so the total contribution is 

oo ^ 

Vphonons = + n;;)- = M/r, (2.13) 

m=2 

where M is the phonon mode number. This sum starts at m = 2 because infinitesimal 
m = 1 modes are equivalent to translations and rotations, and thus are excluded in the 
center of mass frame. The quantization of the string itself produces the universal bosonic 
correction 70 /r first identified by Liischer [29], and so, including the contribution from the 
self-energy, we write the potential as 

PE = V(r) = ^i±J^ + 2nar. (2.14) 
r 

We can now write down our Schrodinger equation (in the p variable) 

Hf = {T + V)f (2.15) 

We recognize this equation as a Sturm-Liouville problem with self-adjoint form [32] 

{d', + Xf-^-^V{pi)}f = 0, (2.17) 

where A/ = ^Ef. 

The traditional Isgur-Paton model looks on the potential ~ 1/r as unphysical for 
r 0. The potential V{r) = {M + jo)/r + 2Trar diverges at small r, and indeed we should 
expect a string model for the flux tube to break down when the radius is on the order of 
the thickness of the flux tube. Thus Isgur and Paton multiplied the offending term by a 
suppression factor 

F{r) = (1 - e"/^'') (2.18) 

with / a free parameter determined by fitting the model to the available lattice spectrum. 
As there are no other adjustable parameters around to influence the spectrum, it was 
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appealing to have something adjustable to minimize x^. However, after studying the 
behavior of the eigenfunctions numerically and analytically, it seems the eigenfunctions 
fall off sufRciently fast as r for the singularity at r — to pose no problem. 

There are other ways to handle the suppression factor. The factor could be applied to 
the whole potential, for instance. Doing so is like supposing that the model ceases to exist 
as r — !■ 0, or that the wavefunction is suppressed by such factor. New evidence coming 
from a study of the dual Ginzburg-Landau flux-ring [24] implies that cr —> as r ^ 0, that 
is, the string (self-energy) disappears as the potential becomes infinite. This approach 
produces some nearly intractable mathematics-some avenues are explored in Chapter [T] 



Chapter 3 

The Isgur-Paton Model in 2D — 
Including Elasticity 

3.1 Canonical variables 

Let's look in detail at the Isgur-Paton model in 2 spatial dimensions. Viewing the flux 
tube as a classical string with 

r = r(0) , [0,27r) (3.1) 

we expand 

r = To + (om cos mcj) + bm sin mcj)) (3.2) 

m=2 

and quantize with respect to ^ to arrive at the one dimensional problem 

T=^-,y = F(r)( —^+2TTar). (3.3) 

47r(Tr r 

Setting F{r) to unity for now, wc examine the kinetic term as we apply the prescription 
for taking a classical system over to quantum mechanics [33]. When P,. we are 

faced with a quandary over what to do with the initial 1 /r coming from the inertial mass. 
One approach would be to multiply through by fi and work with the system 

+ 2,i[V -E]=0, (3.4) 



Q 
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3.1 Canonical variables 



converting the generalized eigenvalue problem to one that is diagonal by means of a 
Cholesky decomposition [34]. Our approach is to find the canonical variables for the 
system. This step actually should be performed at the classical level, before discussion 
of quantum mechanics enters the picture [33], as we want the Hamiltonian to equal the 
total energy, but we can just as readily address the problem after Pr We want to 

transform the independent variable r into some canonical variable p such that the kinetic 
energy becomes proportional to the square of the canonical momentum, ie: 

T{r)^T{p)^-^. (3.5) 

Looking at the troublesome kinetic term, 

p2 _i d_d_ 
T=-^^- — ^r9r here /i = 1. (3.6) 

Symmetric forms of this operator are (upto a constant): 



d fl d 



dr \r dr 



(3.8) 



Evaluating these two operators leads to an ambiguous term in the potential. The first 
operator leads to the equation (for y = 0) 



2 



E = r-'^^r-i (3.9) 
Er^ = ^r-'^ (3.10) 

(3.11) 

(3.12) 

Er=^-~ + "::i, (3.13) 



92 1 




r 2 




d f ,-1 d 


1 


dr \ dr 




1 d^ 


3 d 


r 2 f 

Qj.2 


2 — 

dr 


92 1 d 


3 1 


Q^2 ^ Qj. 


4r2 



3 



As this operator equation acts on 'ip, let ip = fg, then 



f'g + V'g' + fg" - ^ifg + fg') + I^f9 = Erfg (3.i4) 
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3.1 Canonical variables 



Annihilating terms in g\ 



/' _ 1 , 1 /" _ -1 1 



^ = -^/ = ^^V = — 12 (3-15) 



thus 



,311111, , , 

■9[-^-,--^-,--^-,]=Erg (3.16) 



9 = Erg (3.17) 



The second operator gives 



did , , 



dr r dr 
1 92 I Q 

r dr^ r^ dr 



(3.19) 



so 



52 Id 

Er = ^--^. (3.20) 



Q^2 J. Qj, 

Comparing to Eqn 13.131 we note a difference of | ^ . Proceeding, we find 



3 J 

4 r 



^9 = Erg, (3.21) 



and we are still left with a generalized eigenproblem. This ambiguity arises from the choice 
of whether to sandwich the 1 /r factor between the derivatives or vice versa. 

The historical treatment [6] of this term is equivalent to using the first analytical 
form of the kinetic operator, as no contribution to the potential arises from the change of 
variables. Attacking the independent variable directly in Equation 13.61 vields 



^ dr dr 

dr r dr ' 

— pp 
so that conceptually Pr — For our canonical variable p,T , so let 

d dp d 
dr dr dp ' 



(3.22) 



then 



dpd_ 
dr dp 
T 7 
f 2 



(3.23) 



(3.24) 
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so ^ 



^ = r 2 gives 



r^dr = dp => p= -r^, 



(3.25) 



but we want not to change the scale, so we relax our requirement to |^ oc r 2 so that 



p = r^, p{r = 1) = 1. 



(3.26) 



Then 



-1 (iV ^ 

Aira \2j dp^ 
-9 



167rcr dp'^ 



(3.27) 
(3.28) 



and our Hamiltonian becomes 



-9 9^ 
Idna dp' 



+ V{p{r))=E. 



(3.29) 



Absorbing the initial constant into the eigenvalue and potential, 



dp' 



+ \-V 



/ = 



(3.30) 



where p goes from to 00. Dropping the tildes, we follow [35] and solve by first letting 
t = on the interval [0, 1] so that 



m 

dp 



2 dH „ , 
— — 



5p2 



df^ 



2t' 



dt 



(3.31) 

(3.32) 



and then $(t) = tf gives 



df 



(3.33) 
(3.34) 



Our Hamiltonian goes to 



(3.35) 
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Taking t now to be discrete 

tj=jh,h=^^^,je[0,n+l] (3.36) 
and $(f) — > = ^{tj), with the boundary conditions 

$(0) = $(1) = ^ $0 = ^n+i = 0. (3.37) 
The second order approximation 



9*2 ^ 12/i2 



« + 16*i-i - 30$, + 16$,+i - $^+2) (3.38) 



gives us 

Sb^ + ^[V-X]^ = S^ = 0. (3.39) 

The boundary conditions on $' show up on the extreme of the diagonal of the band matrix 
Sb- As we are in 2 dimensions, $„+2 = —^n and $_i = $1. To get the usual form 

Ak = Ax, (3.40) 

use the diagonal matrix K so that 

KSK (A'"i$) = (3.41) 

then 

A - A = '^KSK = '^fSb + V-\ (3.42) 

is the matrix to diagonalize for the eigenstates. {K is the diagonal matrix Kjj = j'^.) 
Orthogonality is given by 



j dpgf ^ (-) I dtt-H^ ^ Yl (3-43) 



If we include the suppression factor F{r) as in [6], we can calculate the spectrum as a 
function of the parameter / in 

F{r) = 1 - e-^'', (3.44) 
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where ^/a = 1 sets the scale in the theory. While the lattice community is moving away 
from the ^/a convention when reporting data, we will continue in ^/a units. The spectrum 
shown in Figure [32] reproduces that of [6,28]. 

3.2 Elasticity 

Looking now at the potential in detail, let's consider the properties a flux tube should 
have. Evidence from the lattice supports the Liischer term for straight flux tubes [30] and 
from type-11 superconductor flux tubes [21] a contribution from the curvature. Nonrela- 
tivistically the curvature term is 



curvature 



7 A 



(3.45) 



For a circular loop, we can think of the Liischer term as a finite additive renormalization to 
the coefficient of elasticity 7 ^ 7 — ^ as we quantize the angular variable (p. Incidentally, 



the string tension also receives a renormalization which we absorb into — > cr- 



renorm ■ 



15r 



M+7 

V = 2jipa + ( -')(1-exp(-f p)) 

, P , 



10 ■ 
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Figure 3.1: The spectrum of the Isgur-Paton model as a function of the parameter 7. For 
comparison, lattice masses are displayed on the right. 
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Our potential, Equation 13. 3[ is now a function of the two parameters 7 and the / from 
Equation 13.441 as well as the variable p. We could now do a comparison of the spectrum 
to lattice data, defining a merit function which we will call (see Equation I3.46P as a 
function of the two parameters / and 7. At first we use the statistical errors as given by 
the lattice data, later adjusting the weightings by manipulating the errors. The behavior 
of the spectrum as a function of 7 at /= 1 is shown in Figure 13.11 The lattice masses 
for SU(2) are given on the right hand side. Its effect is roughly linear, and the signal for 
no curvature term would be 7/^4 = — y| ~ —1. Preliminary work not reported here [36], 
suggests that the suppression F(r) is mostly superfluous, and as we will be adding another 
parameter later to incorporate splitting between the charge conjugate sectors, we might 
as well warm up on a 1 parameter minimization. 

So far the model can produce states with C — + and exhibits parity doubling for states 
with J > 0. Looking at the lattice masses as a function of l/N"^ [1], we see that the 0+ in 
N=2 is a smooth continuation of the 0++ for higher N. This fact will be important later 
as we consider mixing mechanisms to produce the C — — sector. So we select the N=2 
data for our initial comparison, which has only C — + and should test the viability of the 
Isgur-Paton model to explain the lattice spectrum. 

Some ambiguity exists in the assignment of spin to the lattice operators [37]; the 
variational method should select out the state with the greatest overlap with the selected 
symmetry channel and excitation, but convention dictates that the lowest compatible spin 
be assigned. For example, on a cubic lattice the continuous rotational symmetry is broken 
to the discrete 0{2) group-rotations are only defined upto a phase of tt/2, so that spin 
is defined upto mod4. Practically, what is assigned spin could be spin 4; furthermore, 
spin 1 could be confused with spin 3, as it is really |J|mod4 for J = 0, ±l,±2,etc. that 
is determined. Eyeballing the lattice data for N=2 in Table 11.11 and the spectrum in 
Figure we see that the lowest 0+, 2+, audi"*" are in the right ballpark, but the 0~ 
is nowhere near. The lightest state in the Isgur-Paton model with 0~ has M=8 (see 
Appendix IB] for details) and a mass around 15, but there is a state 4~ near the lattice 
mass, suggesting that that operator is coupling to 3—4 and not J=0. To be safe, we will 
only include the well identified states 0+ and 2^ when determining x^- 

For details on why we can use and what it really means, see Appendix[Xl Succinctly, 
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maximum likelihood with uniform priors and normal errors implies minimizing 



(3.46) 



where the sets Mj and Dj are the model's predictions and the data for comparison, and 
the <Tj are the weightings of the data points by the statistical errors. While our statistical 
errors are normal, we expect a systematic discrepancy arising from the breakdown of the 
string picture as r — > 0, so that the smallest state with the tightest statistical weighting, 
the 0+, should be least well represented by the model. We will later deal with this difhculty 
in the next section by introducing a fudge parameter (3, but for now we will stick to using 
the statistical errors and ignore the absolute value of x^, which will be rather large. Our 
problem reduces to one dimensional minimization. Following [38,39], we write the log of 
the posterior 



L = prob{'^\Dk) 
1 2 

= const ^ lyC 



(3.47) 
(3.48) 



Taking a Taylor series about a local minimizer 70, we write 

i = ^(70) + ^1^170 (7-70)' + 



(3.49) 



Stopping at the quadratic term represents a Gaussian distribution centered at 70 with 
variance 

Brute force is a perfectly acceptable method to find the global minimum, so we compute 
X^(7) as shown in the top half of Figure [321 ^^'^ start the numerical minimization routine 
at the minimum of the grid. To get the variance we use the numerical derivative [34] 



/(7 + fe)-/(7-fe) 
2h 



(3.51) 



where the discrete step h is chosen carefully to reduce error. For as in Equation 13.461 



dx^ ^sp Mk- Dk dMk 



97 



97 



(3.52) 
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k 



(3.53) 



Near the minimum 7q we can ignore the second term and write 




1 



cr; 



(3.54) 



3.2.1 Explaing the figures 

The resuhs are shown in hne 1 of Table 15.11 and the bottom half of Figure 13.21 We should 
explain the layout of this figure as it will be repeated throught this thesis. The x-axis is 
given by the total phonon number M, and is related to the spin J. For M = 1 or 2, the spin 
is equivalent. At M — 4 we compare with the lattice state and at M = 5 we compare 
with the lattice states with J = 1. Plotted as a '.' with errorbars are the masses from 
Table [TTTl The ordering of the PC eigenvalues for the lattice states are given in the lower 
left corner. As we are currently in SU(2), the C = — sector is not present. Alongside the 
corresponding C sector are the predictions of the model, plotted as an 'o'. The smallest 

from the grid and the from the fit are presented in the upper left corner. In the 
lower right corner we display the values of the parameters, only 7 in this instance. 

We easily see that the minimum is found at 70 = —.326(141). With only the four 
lightest states at M = and 2, the heavier states are pure predictions of the model. 
The validity of the adiabatic approximation is called into question looking at the excited 
0+. For states with a radius on the order of the width of the flux tube, the adiabatic 
separation of modes is a rather drastic approximation. That so, the agreement which the 
model displays is rather impressive, especially out at M=5 with the J=l from the lattice. 
Looking at M=4, we find a state with J=4^ quite close to the lattice 0^. This feature 
persists throughout the calculations. 
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N = 2 
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Figure 3.2: The posterior and spectrum for N=2. The model's predictions are shown by 
'o', and the corresponding lattice states by '.' with errorbars. The x-axis is labelled by the 
phonon number M. The relationship between M and the spin J is given in Appendix IB| 
Table iBTl See Subsection 13 . 2 . II for more details. 
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3.3 Adjusting the weightings 



Knowing that we have nonnormal systematic effects, we might be interested to what level 
the model is capable of fitting the lattice spectrum. We face the dilemma that the state 
with the greatest weighting in x^, the 0+, is so small that a string picture for the state 
does not really seem applicable. In order to more evenly spread the weighting of the data 
points, we add in quadrature to the statistical error a percentage /3 of the mass of the data 
point Dk, ie 

(3.55) 



The intended effect is to deweight the lightest (and smallest) states so that they do not 
unduly influence the model. So that /3 cx) is not the trivial solution, we must also add 
a regularizing term, thus the quantity to minimize becomes 



a+b=y: 



(3.56) 



To find the extremal value for /3, we take the derivative 



d(3 



{A + B) = 2/3 



1 - 



Mk-Dk 

4 + P'Di 



and hence 



1 



M,, - Dk 



(3.57) 



(3.58) 



4 + P'Dl 

is our equation to solve for (3. While messy, if we fix the Mk to those given by the best 
fit parameters of Equation l3.471 we can solve the equation numerically for /?. Technically 
we should perform a simultaneous minimization over all the parameters [38], but this 
approach will produce a satisfactory qualitative answer. Using the data in Table [5TT] and 
Equation 13. 58i /3 ^ .05, thus we can say the model fits to within 10% (±cr„cw errorbars) 
of the lattice data. 

These new errorbars correct for the systematic discrepancy we expect at small values 
of the radius, where the model is less applicable. Using these new values for the (T^, we 
again perform the minimization of to find 79. Figure [3T3l shows the new posterior and 
the resulting spectrum. The best fit 70 has not changed much, and the value Xfn = 3.39 is 
consistent with 4 data points - 1 parameter. We see that the errorbars on the 0++, while 
still the smallest, are now of the same order of magnitude as the others. Later we will 
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Figure 3.3: The posterior and spectrum for N=2 with adjusted weightings. 



perform parallel fits at /? = and (3 = .05 and compare the results. Essentially the effect 
of /? is to affect the normalization of in 7 space. Over at M=5, the lattice 1+ looks a 
little low (possibly closer to the 3+), and the M=4 is still predicted at the "0~" mass. 



Chapter 4 

Mixing Mechanisms and the 
C = — Sector 



We would like to apply the Isgur-Paton model to N>2, particularly N=3 and the extrapo- 
lation to N=oo. However, these gauge theories posses states with charge conjugation (C) 
eigenvalue of minus [10]. For N>2, the flux tube has an orientation associated with the 
color electric field, giving rise to two orientations L and R which cannot be connected by a 
rotation in two dimensions. Taking ipL,B. to represent these states of opposite orientation 
(see Figure HTTj) . the linear combination ±(-0^ — -0^;) has C = — . These states are degener- 
ate and miique. To go from the degenerate basis [0^, ippi] to [tP+, 0-]; where the subscript 
on tp gives the C eigenvalue, we need a mechanism to transform a into a iJjb. (and vice 
versa) . A transition through an intermediate state with coefficient a immediately suggests 
itself, which we call indirect mixing. We also consider a direct mechanism, whereby at 
small radius the flux tube picture breaks down and a more bag like picture is appropriate, 
allowing a loop of one orientation to flip to the other, with magnitude a. We write the 
Hamiltonian for the enlarged basis as 



Hi 










1 






+ a 









H2 




1 






(4.1) 



where Hi 2 are selected according to the details of the mixing mechanism and operate on 
the relevant subspace of the full Hilbert space of states. 

In the direct mechanism. Hi 2 are simply Hl h which are really each just copies of the 
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(a) ipL (b) ipR 

Figure 4.1: The orientations of the Isgur-Paton flux tube. 



extended Isgur-Paton Hamiltonian Hjp, where 



(4.2) 



thus we can write our problem as 



/ 














Hip 


+ a 


1 






IpL 


V 


Hip 




1 









(4.3) 



For small a, we can regard H' as a perturbation [33], and solve to first order. Then 



\ = Eip± JEjp - {ElEr -a^)=Eip± a. 



(4.4) 



Thus the C=± sectors are spHt evenly away from their mean equal to Eip. The problem is 
that the 0+ ought to be a smooth continuation of the 0++, see [1], not of the average of the 
0++ and . Also, as written, we have not given a any dependence on r. Nonetheless, 
we consider the model viable and calculate its spectrum. 

For the two indirect mechanisms considered, we suppose there exists a state of higher 
mass which interpolates between tpL and ipR. The Hamiltonian is 



; 1 1 1 1 1 the perturbed mass is given by 



( 


Hip 





+ a 


1 






tjj± 







Hi 




1 









(4.5) 



A = ^{Eip + Ei) ± '^^/{Eip-EiY+Aa^. 



(4.6) 
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(a) -ipA (b) i)v 

Figure 4.2: Flux tube states for indirect mixing. 

The state ipi couples to 4'± = '^(''^l ± V'-r) depending on the details of the model. The 
unperturbed eigenstate retains the mass of the Isgur-Paton model, Ejp. 

The vertex mixing mechanism posits a state consisting of the original flux tube plus a 



diameter, see Figure 4.2(b) When N=3, this diameter is a single flux tube meeting the 
circle at two "baryonic vertices" ; for N>3 we suppose the diameter to consist of — 2 flux 
tubes together. The baryonic vertices could contribute to the state's mass directly [8,40], 
but we ignore such contribution for simplicity. To write its Hamiltonian, we include the 
self energy of the additional flux tubes and a contribution from the angular momentum 



Hv = KE+ {2tt + 2{N - 2))ar + —. 



(4.7) 



Here, J is the spin of the state whose mass we are calculating, and I is the classical moment 
of intertia of a ring with N-2 colinear diameters. Specifically, 



{2tt + 2{N -2)/3y 



(4.8) 



The angular momentum also affects the C eigenvalue. Writing 

3 

i^J^Y. e'^"'" ^^6=0 , where - ^/2, (4.9) 

states of even J are invariant under C, while states of odd J transform as C = — . Thus, 
for example, the 0^+ mixes with ipy and is lowered relative to the unmixed , while the 
is given by the unmixed spectrum and the 1^ is lowered. The resulting structure 
might explain the anomaly at J=l in the lattice data: while data for even J shows the 
C = — sector to be heavier than the C = +, for J=l the C = — is nearly degenerate with 
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the C = +. As oo, so does the number of flux tubes m the diameter, and the mass 

should diverge. Therefore we can only do the comparison for the three values of N for 
which we have data. Later we can extrapolate to = oo the values of the parameters if 
we wish. Since there is no analog of V-'y in N=2, the O"'' is at the Ejp mass, and so is not 
a smooth continuation from the mixed O''"^ but of the unperturbed . 

For adjoint mixing, the intermediate state is another Isgur-Paton model flux tube, 



where the flux is in the adjoint representation, Figure 4.2(a) An adjoint flux tube exists 
from N ~ oo down to N=2 [8,41,42], and so the 0+ is a linear combination of i/^/p and 
ipA, and thus is a smooth continuation of the 0+^. The Hamiltonian Ha resembles that 
for Hjp, Equation 14.21 except for the replacement of ap by a a 



(4.10) 



To keep our parameter count down, we accept the hypothesis of Casimir scaling [23,41,43] 
to give a A in terms of crp . Accordingly, the ratio of the string tension in representation j 
to the fundamental string tension equals the ratio of the 2nd Casimir operators [25] , thus 



0_A 

ap 



Ca 
Cp 



2N^ 

Af2 - 1 



(4.11) 



is the value we take for each N. 

At this stage, we are in order to compute the log of the posterior (Appendix [A| for 
each model 



L = lnpro6(7, a\{Dk}, {cfc}, /) 



(4.12) 
(4.13) 



in the Gaussian approximation. Now that we have two parameters, we expand 



L = L{jo,ao) + - 



d'^L d'^L d'^L 

^l7o,ao(7 - 7o)^ + -^ko.aol" - "o)^ + 2^:^l7o,ao(7 " 7o)(a - ao) 

(4.14) 



Isolating the quadratic part, we write the curvature matrix 



[C] = 



d'jda. 



d'jda do 



(4.15) 
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whence the covariance matrix is 



(4.16) 



As 



then 



and 



^7 2 d'y ' 



jk 



1 a^x^ 



djdk 



2 ajafc 



where j, A; e {7, a}, 



(4.17) 



(4.18) 



(4.19) 



dj dk djdk 

Near the minimum of x^, where the expansion is happening, we can ignore the nonUnear 
term and write 



1 



dMj \ dMj dMi 
(97 1 c?7 da 



dMj dlvh ( dMj \ 2 
da \ da J 



(4.20) 



which is evaluated by 



92/ 



d^da 



[/(7 + /I, a) - /(7, a) - /(7, a + /i) + /(7, a)]/^^ 



(4.21) 



for some appropriately chosen h [34]. To find (70,0:0), we again apply brute force to a 
grid in parameter space to map 'X^{'y,ct) and begin a simplex [34] minimization from the 
minimum on the grid. 



Chapter 5 



Analysis I 

5.1 Fitting the data 

Let's begin by using the statistical errors, i.e. (3 — 0: these fits are driven mostly by the 
lightest state, the 0++. For direct mixing, as the flux tube loses its orientation for N=2, 
we should start by looking at Figure from the preceding chapter. By using only the 0+ 
and the 2+ states when fitting the parameter, states at higher M are predictive. We see 
that the 1^ agrees within la, but the 1+ is as close to a state with M=3 as it is to M=:5. 
The question of lattice spin assignments is further called into question by the 0". Here, 
the model predicts a state with M=J=4 quite close to the lattice 0~. The lightest state 
with 0^ quantum numbers in the Isgur-Paton model has M=8 and has mass around 16 
in our units, much too heavy to match the lattice estimate. The first excitation of 0^ is 
included in the fit. The model's prediction is too high by about 4cr. We should remember 
that the model's adiabatic approximation separating the radial from the phonon modes 
has little justification and will contribute to the systematic error. 

When N=3 we encounter our first 2 dimensional parameter estimation problem. The 
posterior density function (pdf) is shown in Figure 15.11 Despite the coarse graining of 
parameter space, we see a clearly defined maximum in the pdf. The best is at 7 = 
0.54, a — 1.89. Examining the best fit spectrum, we see a qualitatively good fit. Except 
for the 2^ , everything used in the fit is within 2(t of the lattice values. As the mixing 
is just a constant, no properties of the states can infiuence the magnitude of the mixing, 
which might account for the discrepancies at higher M. These states are larger and thus 
less likely to collapse to the ball pictured for this mechanism. The discrepancy in the 
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C = — sector continues for the M=4 and 5 states. Here we see good predictive agreement 
for the C — + sector (ahowing for the spin assignment at M=4), but the C ~ — sector 
does not agree at higher M. If we are concerned with an 0/4 ambiguity, then we also have 
to consider a 1/3 ambiguity. In that case, there is a J=M=3*~ at the right focation in 
the spectrum, but then the 3++ should couple to what is assigned to J=l. 

What should we make of the mess in the heavier part of the spectrum? First, we should 
stress that these states are not included when determining f^nd hence the parameters 
- they are free predictions of the model, which is what raises the question of ambiguities 
in the conventional assignment of spin to certain lattice operators. A preliminary explo- 
ration [44] showed that for the ++ sector the J=0 should be distinguished from the J=4. 
But states with more complicated quantum numbers and higher spin require operators of 
more complicated construction, and the situation regarding which states a (finite) opera- 
tor couples to is much less clear [37]. It is possible that the 1*+ operator is good while the 
is being misled by the J=3, but the issue requires further investigation. Exploratory 
calculations on a new technique of constructing lattice operators for arbitrary spin are 
discussed in Chapter [HI 

For N=4 and 5, the qualitative features of the spectrum do not change dramatically. 
The low lying states reproduce the lattice spectrum reasonably well, and the higher states 
continue to follow the spin ambiguity predictions. Turning to the posterior. Figure 15. 3[ 
at N=:5 we start seeing evidence of a second peak near 7 = 0. This peak is not included 
when we take the local quadratic approximation to compute the errors on 7 and a, but 
its presence should be noted. The main peak has also shifted closer to the 7 = axis. 

The lattice spectrum has been extrapolated to N=oo using a correction linear in 
(1/A^^) [1]. Repeating our fitting procedure with these lattice values, we get the posterior 
and spectrum of Figures l5.4l Now the second peak has swamped the first to dominate the 
posterior, though there is a little of the first peak left. The best fit gives 7 — —.042(190) 
with an error consistent with 7 = 0. Were 7 to equal zero would require a delicate balance 
of the elasticity and the Liischer correction to cancel. We will look at the behavior of 7 
as a function of (1/A^^) later in the analysis and compare its extrapolation to N=oo with 
7 fit to the extrapolated spectrum. 

Looking now at the adjoint mixing mechanism, we can start with the N=2, as there 
still exists an adjoint loop mixing with the fundamental loop, which is plotted in the 
figure. Our prior nearly misses the peak. Figure 15.51 but as its location moves into the 
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prior for later N, we should feel confident we are seeing the global minimum. The low spin 
spectrum fits quite well - all states are within 1.5a. The higher spin states also display 
good agreement. 

At N=3, the model starts showing signs of difficulty in fitting the lattice data. Including 
both radial excitations and the C — — sector contributes many more states to the fit. Still, 
we see a well defined peak in space and can determine 70 = .576 and ao — 3.56. Looking 
at its best fit spectrum, the model is in general qualitative agreement with the lattice data. 
The most obvious discrepancy is the 2^^. At higher spins, the model consistently is too 
heavy. In the C = — sector, the oddly low 1=*=^ does compare well with the 3^^. The 
magnitude of the radial splitting seems a little off as well, again perhaps because of the 
adiabatic approximation. As — > 00, nothing dramatic happens. The values of the best 
fit parameters. Table [HTTl do change somewhat. The best fit 7 are still positive, and a ^ 3. 
Except for the lowest J=0 state in either C sector, the model's predictions are too heavy, 
especially at higher spin. 

For the vertex model we have only three cases to look at. From N—3 to N=5, we 
find a positive 7 and a mixing a ~ 2, Figures [5. 101 through 15.1^ For J^O and 2, the 
model generally reproduces the features of the lattice spectrum. The robustness of the 
Isgur-Paton model is shown in the consistency of results across mixing mechanisms. At 
higher spin, things are a mess, except for 4 ^. That is, if we persist in comparing the 
lattice J=l state to M=5. Allowing for more spin ambiguity, we might compare these 
states to M=3. In that case, the 3^+ would be in good agreement, while the 3*~ might 
also be as well; the lowest state 3^^ in the model is mostly ^pv, while the lattice operator 
is not constructed to couple to such a state [45]. At this point, such remarks are only 
conjecture - a much more thorough technique needs to be developed, as in Chapter m As 
N 00, the features remain the same. One state that does not fit well is the 0^ , which 
the lattice claims is well below the model's prediction for the 4*^. 
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5.1 Fitting the data 



N = 3 





Figure 5.1: The direct mixing posterior and spectrum for N=3. The model's predictions 
are shown by 'o', and the corresponding lattice states by '.' with errorbars. The x-axis is 
labelled by the phonon number M. The relationship between M and the spin J is given in 
Appendix [Bl Table IbTTI See Subsection 13.2. II for more details. 
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N = 4 
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Figure 5.2: The direct mixing posterior and spectrum for N=4. 
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5.1 Fitting the data 



N = 5 
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Figure 5.3: The direct mixing posterior and spectrum for N=5. 
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5.1 Fitting the data 
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Figure 5.4: The direct mixing posterior and spectrum for N=(X). 
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5.1 Fitting the data 



N = 2 





Figure 5.5: The adjoint mixing posterior and spectrum for N=2. 
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5.1 Fitting the data 



N = 3 
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Figure 5.6: The adjoint mixing posterior and spectrum for N=3. 
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N = 4 
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Figure 5.7: The adjoint mixing posterior and spectrum for N=4. 
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N = 5 
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Figure 5.8: The adjoint mixing posterior and spectrum for N=5. 
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Figure 5.9: The adjoint mixing posterior and spectrum for N=oo. 
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5.1 Fitting the data 



N = 3 
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Figure 5.10: The vertex mixing posterior and spectrum for N=3. 



39 



5.1 Fitting the data 



N = 4 
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Figure 5.11: The vertex mixing posterior and spectrum for N=4. 
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N = 5 
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Figure 5.12: The vertex mixing posterior and spectrum for N=5. 
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5.1 Fitting the data 



N direct 


0++ 0— 


2±+ 2±- 


3±+ 3±- 


4±+ 4±- 


1±+ 1±- 


2 rrtjpc 

lll'jj.r 


4.618 
7.482 


7.946 
9.9G3 


9.324 
11.09 


10.56 
12.13 


11.67 
13.11 


3 


4.278 6.171 
6.707 8.60 


7.261 9.154 
9.051 10.94 


8.509 10.4 
10.11 12.00 


9.638 11.53 
11.09 12.99 


10.67 12.57 
12.02 13.91 


4 


4.146 6.072 

6.601 8.526 


7.151 9.076 

8.955 10.88 


8.407 10.33 

10.02 11.94 


9.542 11.47 
11.01 12.93 


10.58 12.51 

11.93 13.86 


5 


4.062 5.715 
6.612 8.265 


7.145 8.798 
8.998 10.65 


8.431 10.08 
10.07 11.73 


9.589 11.24 
11.08 12.73 


10.65 12.3 
12.02 13.67 


oo 


4.009 5.144 
6.279 7.853 


7.219 8.354 
9.155 10.29 


8.552 9.687 
10.26 11.39 


9.747 10.88 
11.28 12.42 


10.83 11.97 
12.24 13.38 



Table 5.1: The spectrum for direct mixing at various N. For each N, the first mass is the 
lowest state in that channel, and the second mass is the radial excitation. 



N adjoint 


0++ 0~ 


2±+ 2±- 


3±+ 3±- 


4±+ 4±- 


1±+ 1±- 


2 rrijpc 
Tn*jpc 


4.717 
7.04 


7.643 
9.432 


8.888 
10.51 


10.02 
11.51 


11.07 
12.46 


3 


4.217 6.225 
6.801 8.64 


7.417 9.196 
9.329 10.98 


8.754 10.44 
10.46 12.03 


9.961 11.57 
11.51 13.02 


11.06 12.6 
12.49 13.94 


4 


4.083 6.134 
6.697 8.572 


7.316 9.124 

9.244 10.92 


8.664 10.38 
10.38 11.98 


9.879 11.51 
11.43 12.97 


10.99 12.54 

12.42 13.89 


5 


3.999 5.813 
6.687 8.336 


7.322 8.874 
9.283 10.71 


8.698 10.15 
10.43 11.79 


9.934 11.3 
11.5 12.79 


11.06 12.36 
12.5 13.72 


oo 


3.919 5.78 
6.621 8.312 


7.258 8.849 
9.226 10.69 


8.639 10.13 
10.38 11.77 


9.879 11.28 
11.45 12.77 


11.01 12.34 
12.45 13.71 



Table 5.2: The spectrum for adjoint mixing at various N. For each N, the first mass is the 
lowest state in that channel, and the second mass is the radial excitation. 



N vertex 


0++ 0— 


2±+ 2±- 


3±+ 3±- 


4±+ 4±- 


1±+ 1±- 


3 rrijpc 
ni*jpc 


4.4 6.036 
6.742 8.5 


6.839 9.049 
8.748 10.86 


10.31 7.975 
11.92 9.719 


9.007 11.45 
10.62 12.91 


12.49 6.765 
13.84 8.849 


4 


4.203 6.014 

6.735 8.484 


6.878 9.031 

8.905 10.84 


10.29 8.113 

11.91 9.961 


9.238 11.43 

10.94 12.9 


12.47 7.448 
13.83 9.649 


5 


4.074 5.606 
6.762 8.185 


7.006 8.713 
9.102 10.58 


10.01 8.331 
11.66 10.22 


9.531 11.17 
11.26 12.67 


12.24 7.947 
13.61 10.28 



Table 5.3: The spectrum for vertex mixing at various N. For each N, the first mass is the 
lowest state in that channel, and the second mass is the radial excitation. 
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5.2 Adjusting the weightings 



5.2 Adjusting the weightings 

How do things change when f3 = .05? Actually, not by very much. By adding ±5% to the 
error, we de-weight the states with the tightest statistical errors, which are the low lying 
states of small radius on the order of the thickness of the flux tube, allowing the radial 
excitations to take on more of a role in determining the best fit parameters. We are also 
trying to quantify what we can expect of a nonrelativistic model with a drastic adiabatic 
approximation when compared to our best estimates of the pure glue spectrum from the 
lattice. 

Again considering direct mixing first, we look at the log of the posterior as a function 
in parameter space, Figure 15.131 We note a broader, smoother peak located at roughly 
the same location (70, ao). Consistency in the value of the best fit parameters implies 
that we have not done anything too drastic by manipulating the weightings. A broader 
posterior actually implies a more conservative estimate of the parameters, in the sense of 
having a larger errorbar (see [38,39]). The spectrum, for N=3, Figure [5?T^ shows good 
agreement across the lower states - except for the and the 2 ^, all the predictions 
are within the expanded errorbars. The higher states are again consistent with ambiguous 
spin assignments at M=:4 and 5. As iV ^ 00, Figures [5.141 through 15.161 the features 
of the model are essentially the same. The prediction for the is too low, which 
could also be explained if the tighter errorbar of the 0++ were pulling the whole spectrum 
towards the state which it should least fit. At N=oo, where we are fitting parameters to 
an extrapolation of the lattice data, our fit for 70 = —.046(295) is negative, but within 
its error is consistent with zero. We will more closely examine the behavior of 7 and a as 
functions of N in the next chapter. 

For the adjoint mixing mechanism at N=2, we see a very elongated posterior (and a 
miniscule value for x^)- The peak is still in about the same location, however, which is 
good. The length of the peak tells us that a wide range of parameter values are consistent 
with the lattice data. When N=3, our peak takes on a more traditional shape at 70 = 
.44(36) ; ao — 3.33(60). Looking at the spectrum. Figure [5.171 we see good agreement 
for J=0 and 2 and the usual mess for the higher spins. As 00, the 1=*=" states 

are nowhere close, suggesting either they are spin 3 or that the model does not have the 
necessary structure to explain the data. 

As before, we have only three cases for vertex mixing. For each N, the peak is well 
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5.2 Adjusting the weightings 



defined, giving parameter values in Table 16.21 The spectrum again demonstrates the 
mechanism's complex structure. At M=5, the direction of the C splitting is in the right 
direction to follow the lattice assignments, but far too large in magnitude. These ambigu- 
ities highlight the need for better determination of the spin assignments. 



N direct 


0++ 0— 


2±+ 2±- 


3±+ 3±- 


4±+ 4±- 


1±+ 1±- 


2 nijpc 
Tn*jpc 


4.519 
7.413 


7.869 
9.902 


9.255 
11.03 


10.49 
12.08 


11.61 
13.06 


3 


4.206 5.948 
6.693 8.435 


7.237 8.98 
9.058 10.8 


8.504 10.25 
10.12 11.87 


9.647 11.39 
11.12 12.86 


10.69 12.43 
12.05 13.79 


4 


4.005 5.758 
6.543 8.296 


7.079 8.832 
8.926 10.68 


8.361 10.11 
10 11.75 


9.516 11.27 
11 12.76 


10.57 12.32 
11.94 13.69 


5 


3.959 5.535 
6.558 8.134 


7.082 8.658 
8.961 10.54 


8.383 9.958 
10.05 11.62 


9.552 11.13 
11.06 12.63 


10.62 12.19 
12 13.58 


oo 


3.925 5.137 
6.35 7.849 


7.137 8.349 
9.073 10.29 


8.47 9.682 
10.17 11.39 


9.666 10.88 
11.2 12.41 


10.75 11.97 
12.16 13.37 



Table 5.4: The spectrum for direct mixing with adjusted weightings at various N. For each 
N, the first mass is the lowest state in that channel, and the second mass is the radial 
excitation. 



N adjoint 


0++ 0— 


2±+ 2±- 


3±+ 3±- 


4±+ 4±- 


1±+ ti- 


2 mjpc 
m J PC 


4.714 
7.075 


7.68 
9.486 


8.938 
10.57 


10.08 
11.58 


ll. 14 

12.53 


3 


4.086 5.987 
6.732 8.463 


7.357 9.01 
9.298 10.83 


8.717 10.27 
10.44 11.89 


9.94 11.41 
11.5 12.88 


11.06 12.46 
12.49 13.81 


4 


3.869 5.82 
6.568 8.341 


7.2 8.88 
9.169 10.72 


8.58 10.16 
10.32 11.79 


9.82 11.31 
11.39 12.79 


10.95 12.36 
12.39 13.73 


5 


3.826 5.577 
6.573 8.164 


7.227 8.691 
9.218 10.56 


8.628 9.987 
10.38 11.65 


9.882 11.15 
11.46 12.66 


11.02 12.22 
12.47 13.6 


oo 


3.733 5.422 
6.511 8.052 


7.186 8.57 
9.193 10.46 


8.602 9.879 
10.37 11.55 


9.867 11.06 
11.46 12.57 


11.01 12.13 
12.47 13.52 



Table 5.5: The spectrum for adjoint mixing with adjusted weightings at various N. For 
each N, the first mass is the lowest state in that channel, and the second mass is the radial 
excitation. 



N vertex 


0++ 0— 


2±+ 2±- 


3±+ 3±- 


4±+ 4±- 


1±+ 1±- 


3 mjpc 
m jpc 


4.456 5.825 
6.75 8.345 


6.95 8.884 
8.869 10.72 


10.16 8.106 
11.79 9.863 


9.15 11.31 
10.77 12.79 


12.36 6.678 
13.73 8.85 


4 


4.111 5.701 
6.682 8.254 


6.879 8.787 
8.931 10.64 


10.07 8.147 
11.72 10.01 


9.296 11.23 
11.01 12.72 


12.29 7.28 
13.66 9.609 


5 


3.992 5.469 
6.702 8.086 


6.962 8.607 
9.071 10.49 


9.912 8.3 
11.58 10.2 


9.51 11.09 
11.24 12.6 


12.16 7.848 
13.54 10.24 



Table 5.6: The spectrum for vertex mixing with adjusted weightings at various N. For 
each N, the first mass is the lowest state in that channel, and the second mass is the radial 
excitation. 
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N = 3 
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Figure 5.13: The direct mixing posterior and spectrum for N=3. 
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N = 4 
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Figure 5.14: The direct mixing posterior and spectrum for N=4. 
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N = 5 
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Figure 5.15: The direct mixing posterior and spectrum for N=5. 
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Figure 5.16: The direct mixing posterior and spectrum for N=oo. 
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5.2 Adjusting the weightings 



N = 2 
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Figure 5.17: The adjoint mixing posterior and spectrum for N— 2. 
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5.2 Adjusting the weightings 



N = 3 
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Figure 5.18: The adjoint mixing posterior and spectrum for N=3. 
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5.2 Adjusting the weightings 



N = 4 
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Figure 5.19: The adjoint mixing posterior and spectrum for N=4. 
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N = 5 
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Figure 5.20: The adjoint mixing posterior and spectrum for N— 5. 
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Figure 5.21: The adjoint mixing posterior and spectrum for N=oo. 
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N = 3 
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Figure 5.22: The vertex mixing posterior and spectrum for N=3. 
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N = 4 
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Figure 5.23: The vertex mixing posterior and spectrum for N=4. 
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N = 5 




16 



10 









1 

12.1 






1 


1 


1 

X 






11.4 






X 


X 


O 

X 










ir 


X 


O 


f II 


Il°- 




1° 

X 


ir 

Hx 


ir 




r 

o 




o - 




































++ -+ — +- 

1 






1 


0.139 

1 


2.46 

1 





-1 1 2 3 4 5 6 

M = i; m (n* + n ) 

^ m m' 



Figure 5.24: The vertex mixing posterior and spectrum for N=5. 



Chapter 6 

Analysis of the Parameters and 
the Merit of the Models 



Following the observations of Teper [1], the value of physical quantities derived from 
lattice calculations should demonstrate a dependence linear in The most obvious 

example is the masses themselves - the fits shown with N — oo are actually fit to the 
lattice spectrum as extrapolated to = oo in [1]. Having lattice data at several values 
of N is what allows the extrapolation to be made. As we have collected values for the 
parameters at several N, we can do the same extrapolation to investigate how well the 
parameters follow the (1/iV^) dependence. We take the values at N=3, 4, and 5 when 
performing the fit for two reasons: for consistency across the models and so that the 
values at N=2 are predictive as well as those at = oo. In Table I^TTl the last row for 
each model is the parameter values as extrapolated by (l/iV^)|jv=oo- 

To perform the extrapolation, we find the best m and c in a least squares sense when 
minimizing 

2^ -2 ' (6-1) 

Af=3 



which can be evaluated directly from the data [38,39]. Figure 6.1(a) shows the results 
for 7 at /3 = 0. Across the range of N, a (1/iV^) dependence is quite believable. The 
adjoint model in particular shows a good dependence from N=:2 down to = oo. The 
direct model docs not agree well at the endpoints, but the estimates of "fN=oo from all 



three models agree within errors at 7 ^ .25. For f3 = .05, Figure 6.1(b) much of the same 



behavior is observed. The value of 700 drops slightly to be ~ -1(3), which is consistent 
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with zero within errors. As we are seeing {l/N"^) dependence, we should beheve that we 
are measuring a physical parameter. 

For the mixing parameter a, we see more variety among the models. At /? = 0, third 



column of Table [01 and Figure 6.2(a) we see dependence, but the two indirect 

mixing models have slopes of opposite sign. However, the two models reach a common 
value at iV = oo of 7o ^ 3. That these two alternatives reach agreement at the N oo 
limit suggests that this form of mixing mechanism is physically viable. If these more 
exotic states exist in the pure gauge theory, we should construct a more generalized model 
a la [46], with a Hamiltonian 



H = 



H+ 





ai 


as 







a2 




ai 




Hv 




as 


a4 


0:5 


Ha 



(6.2) 



but that work has yet to be done. The direct mixing mechanism gives a reasonable 
(1/A^^) fit, with an intercept lower than for the indirect mixing mechanisms. As indirect 
mixing is second order in a, and the direct mixing is first order, it seems plausible that 
^direct ~ -y/aindircct at iV = oo. When /3 = .05, thc magnitude of the parameters decreases 



slightly, but the same behavior is observed (see Figure 6.2(b)). 
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N 


7 


(0-7) 


a 






direct 












2 


-.33 


(14) 








3 


0.54 


(20) 


1.89 


(27) 


(23) 


4 


0.49 


(22) 


1.93 


(29) 


(24) 


5 


0.28 


(26) 


1.65 


(34) 


(29) 


oo 


-.042 


(190) 


1.14 


(26) 


(20) 


(ViV2)|oo 


0.23 


(13) 


1.67 


(22) 




adjoint 












2 


1.20 


(47) 


4.55 


(72) 


(58) 


3 


0.58 


(20) 


3.56 


(33) 


(26) 


4 


0.52 


(22) 


3.57 


(36) 


(28) 


5 


0.33 


(26) 


3.15 


(44) 


(34) 


oo 


0.32 


(31) 


3.19 


(52) 


(40) 


(l/A^')|oo 


0.28 


(13) 


3.15 


(36) 




vertex 












3 


0.46 


(20) 


1.88 


(29) 


(24) 


4 


0.45 


(22) 


2.52 


(31) 


(25) 


5 


0.22 


(26) 


2.54 


(39) 


(31) 




0.20 


(13) 


3.07 


(27) 





Table 6.1: Parameter values for /3 = 0. 



N 


7 


(0-7) 


a 


(o-a) 




direct 












2 


-.38 


(30) 








3 


0.41 


(36) 


1.74 


(49) 


(40) 


4 


0.30 


(36) 


1.75 


(49) 


(40) 


5 


0.18 


(37) 


1.58 


(49) 


(40) 


00 


-.046 


(295) 


1.21 


(41) 


(29) 




0.079 


(302) 


1.56 


(54) 




adjoint 












2 


1.09 


(93) 


4.30 


(1.39) 


(1.14) 


3 


0.44 


(36) 


3.33 


(60) 


(45) 


4 


0.34 


(36) 


3.31 


(60) 


(45) 


5 


0.20 


(37) 


2.96 


(62) 


(46) 


00 


0.11 


(39) 


2.82 


(66) 


(50) 


(l/A^')|co 


0.10 


(30) 


2.89 


(84) 




vertex 












3 


0.34 


(35) 


1.59 


(53) 


(41) 


4 


0.27 


(36) 


2.24 


(53) 


(42) 


5 


0.14 


(36) 


2.46 


(56) 


(43) 




0.066 


(302) 


2.98 


(68) 





Table 6.2: Parameter values for /3 = .05. 
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How do we evaluate the merit of our models? While rules of thumb abound, the 
rules of inferential statistics [38,39] highlight the difhculties. Writing the posterior for the 
probability of model M given a set of data D with errors a as 

MMP..)^ '''^';!"'^"' . (6.3) 

the critical factor p{D,a) (named the "evidence") is incalculable. What is the value of 
"the probability that this particular data set was drawn from the myriad possibilities of 
the universe"? When certain simplifying assumptions are met, the likelihood should follow 
the statistical noise, ^ but with few data points and/or nonnormal (systematic) 
errors, the magnitude of a "believable" becomes hard to calculate. Of course, the right 
answer to the question "Is this a good model?" is the question "Compared to what?" 
Naturally, without alternatives we should believe the only model we have! This ambiguity 
disappears when we compare the relative probabilities of two (or more) models 

piMi\D,a) ^ p{D,cT\M,) ^ p(Mi) 

p{M2\D,a) p{D,a\M2) piAh)' ^ ' ' 

as the evidence cancels between the numerator and the denominator. The second factor 
is the ratio of priors, quantifying any prejudice on the models we might wish to incorpo- 
rate and usually set equal to unity. So, to compare two models we should compare the 
(exponential) of their x^ . In detail, there are also factors arising from the priors (which 
cancel in our case, being equivalent) as well as the breadth of the posterior as embodied 



in <T^ and CTq 



p{Mi\D,a) a^iaai 

— ■ exp 



-^ixl ~ xl) 



(6.5) 



p{M2\D,a) a~f2<Ja2 

If we want the factor by which one model is more likely than another, we can take the 
exponential of the difference in x^ ■ For our purposes, we can look directly at x^, our merit 
fimction. Table 

The first column is for /? = 0. The jump from N=2 to N=3 for adjoint and direct 
mixing results from the addition of more states to the fit. The reduction in x^ for N=cx) is 
because of the larger errorbars on the extrapolated lattice data. Looking at the physically 
interesting case of N=3, the direct mixing comes in with x^ — 55.9, while the adjoint has 
X^ — 78.2 and vertex mixing has x^ ~ 85.3. Thus we must conclude that the data support 
direct mixing over indirect mixing by a relative probability of ~ ea;p(— 55.9+80) « 3x 10^*^. 
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N 


/3 = 


/? = .05 


direct 






2 


49.1 


3.39 


3 


55.9 


6.79 


4 


51.5 


9.9 


5 


43.6 


9.56 


oo 


24.2 


8.42 


adjoint 






2 


8.43 


.777 


3 


78.2 


8.87 


4 


79.1 


14.1 


5 


72.6 


13.4 


oo 


47.7 


12 


vertex 






3 


85.3 


8.04 


4 


53.1 


9.72 


5 


55.5 


11.4 



Table 6.3: for th^ Isgur-Paton model 

Looking at the N=oo fits, while the magnitude of the changes, the relative difference 
^ —24.2 + 50 is about the same. Comparing the indirect mechanisms, we note that for 
N=3, 4, and 5 the adjoint's is between 70 and 80, while the vertex method's drops 
to 50 for N=4 and 5. At N=4, the vertex and the direct mechanisms have different 
by only 2. 

After we have adjusted the weightings by /3 = .05, we get x^ ^-t ^ much more familiar 
magnitude. Most of the fits have 12 data points - some have 11. With two parameters, 
that leaves ~ 10 degrees of freedom. We can now apply our rule of thumb x^ ^ d.o.f . to 
the question "Do we believe this model can explain the data at a 10% level?" Looking 
again at Table [^751 and ignoring N=2, the lowest x^ is for direct mixing at N=3, equal 
to 6.79. As the effect of finite /3 is to de- weight the 0++, we should conclude that the 
Isgur-Paton model can explain both the C = + and C = - sectors of the lattice data 
to a reasonable level of agreement. For higher N, the indirect mixings do not fit as well 
as the direct mechanism, and the vertex mechanism is slightly better than the adjoint. 
Even at the 10% level, the indirect mechanisms cannot quite be said to be "believable", 
though they do capture certain qualitative features of the spectrum. Trying to capture 
the complicated dynamics in a simple nonrelativistic semi-classical Hamiltonian might be 
an impossible task. 

We might try to improve the model by making it relativistically correct by including the 
chromo-magnetic contributions. Another avenue is to explore what happens as r — s- and 
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the string picture breaks down. Shedding hght here is a recent study of the dual Ginzburg- 
Landau flux ring [24]. The dGL picture might be applicable to the QCD vacuum, if the 
conditions of Abelian dominance and monopole condensation hold. Defining their effective 
string tension as the total rest energy of a solution of the differential equations divided by 
its circumference, those authors find 

feff — > as r ^ 0. (6.6) 

The traditional Isgur-Paton model has the constant a = IVr, thus these two pictures have 
radically different viewpoints on the properties of a flux tube at small radius. 



Chapter 7 



Analysis II: Including a cr^fj 

One problem for the Isgur-Paton model is understanding what happens as r ^ 0. The 
string picture of a flux tube is only appropriate on scales where the radius is much larger 
than the thickness of the flux tube ^ With the radius ^ -k=, the flux tube would look 
more like an annulus in two spatial dimensions or a torus in three dimensions, and the M /r 
potential for the phonon contribution seems unreasonable. Indeed, the fluctuations might 
not even be able to be described by a phonon potential. A detailed theory of membrane 
dynamics might address the problem, but such is beyond the scope of this work. 

What if as r — > 0, the string itself "dissipated"? A study of the dual Ginzburg-Landau 
theory of type-II superconductors [24] determined a function <Tcs{f')^ Figure [7TT1 defined 
to be the total energy of a flux ring solution of radius r divided by the circumference. 
The total energy was calculated from the microscopic differential equations for the fields. 
The non-abelian gauge theory's flux tube is often thought to arise from a dual Meissner 
effect for a type-II superconducting picture of the QCD vacuum [19,20]. As a is defined 
from the interquark potential as r oo, its behavior as r ^ is uncertain, and we might 
apply fj — !■ (Tcff(r) to the Isgur-Paton model. Transcribing the function from Figure FTTl 
and converting scales, we can fit the function with the form 

^cft - 1 - e-f-' (7.1) 

for f=1.7177, best in a least squares sense. Knowing the functional form of <Toff(?'), we can 
investigate its effect on the Isgur-Paton model. 

When we quantized the model earlier, the difficult term was the kinetic energy, so that 

fi4 
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^,"0.5 



= o~u G D ' O o o 



o from hep-ph/9902441 
1-exp(-1.7177 * r) 



(1-exp(-2.696* r)r 



1.5 2 



2.5 



3.5 



Figure 7.1: CTeff vs r 



is where we will begin. Abusing notation slightly, 



p2 

rp r 

~ 2/i 
-1 

4-77 



ar 



(7.2) 
(7.3) 



where a is now included in the r dependence of the kinetic energy. If a were still a constant, 
we would get (as before) 

-1 9 

—-— (7 4) 

47raconst4 V' ' 



since dp = r^dr, p = ^r^ . Now we must include cjefi in our differential equation relating 
r and p. For some canonical variable p, 
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so 

dp = dra^gr^ (7.6) 

with boundary conditions 

p(0) = 0, p(l) = 1. (7.7) 
Substituting our function for acs, we find 

p = J drr^l-e-f'')K (7.8) 

This integral cannot be put into cfosed form, yet we need the antiderivative explicitly 
to substitute into the potential V{r{p)). While it would be nice to include explicitly the 
function from Eauation l7.1i any similar function which captures the salient features should 
suffice to test its influence on the spectrum. We need a function which goes to zero as 
r —^ and which approaches unity asymptotically, and which is integrable with respect 
to the measure r^^dr. After looking at several functional forms, the closest one meeting 
these criteria is 

aeff= (l-e-^''-)', (7.9) 
where f is now given by 2.696. Performing the integral (and dropping the prime) 

j drri (1 _ e-f^) = ^ (^e'f^ + ^/r) - ^crf /f^, (7.10) 

thus 

p^-ri + ^e-f^'-^eiiJ]^. (7.11) 
Applying the condition p{r = 1) = 1, 

p(l) = c = 0.495544, (7.12) 

so 

1. / 2 „ \ 

C fc 

We cannot find the functional inverse explicitly, but as ultimately we are working with a 



p^lp^y-[ e-f^ + tt/H - -^erf V/r. (7.13) 
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900 1000 



Figure 7.2: r vs p 



discrete set {pj} related to {tj}, we can solve the equation numerically to give rj 



n+1 



(7.14) 



which we can put into the potential 



V = V{r{p,)). 



(7.15) 



The solution of Equation (|7.13[) is shown in Figure 17.21 Returning to the kinetic energy, 



(7.16) 



so our Hamiltonian goes to 



T + V = E^ 
+ 4:nc^V{r{p)) = Attc^E = A, 



(7.17) 
(7.18) 



which we transform as before to the matrix equation 



(7.19) 



and solve. 

So, how does this model perform? Looking at Table [TTTl we see the values of for 
various N at /? = 0. For the direct mechanism, decreases from 117 at N=3 down to 
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N 


2 3 4 5 oo 


direct 


117 79.3 36.5 15.3 


adjoint 


13.2 135 101 81.1 26.5 


vertex 


257 130 96.4 



Table 7.1: for (Teff 

15.3 for N=cx). This N dependence of seems a little odd. For small N, is too large 
to believe the model including the effective string tension, but as — s- oo, x^ drops to a 
value which is nearly believable, even with the tight errors. As the details of the conjugate 
splitting might depend on N, we are seeing evidence that the Isgur-Paton model is most 
applicable to the — > oo limit. 

Looking now at the indirect mechanisms, we see the same type of behavior or . For 
the adjoint mechanism, x^ ranges from 135 for N=3 down to 26.5 for N=oo. Apparently, 
having Ccff — > works reasonably well for — > oo, but for small N including such 
behavior is detrimental to the model. For vertex mixing, x^ simply remains unbelievably 
large, > 100. As this mechanism has an extra flux tube for its diameter, imposing (Toff — > 
is a major complication to the model and is not supported by the evidence. 

Comparing the values of x^ at /? = for both a = \ and (Teff — > 0, we can say that the 
lattice data best supports the direct mixing mechanism in the limit as TV — > oo. Were we to 
believe in using the incomplete gamma function [34] to assess this model's goodness-of-fit, 
we would find 



Q 



r(f) 



(7.20) 
(7.21) 



Turning now to the parameter values for the direct and adjoint mechanisms. Table [721 
we find some interesting behavior. The (l/A'^^) dependence is not readily apparent for a 
and is dubious for 7. Overall, the results seem to split into two sectors across mechanisms, 
one for N=:3 and 4, the other for N=5 and 00. The smaller N return a 7 which is 
positive ~ .6 or .7, while for higher N a negative 7 is supported ^ -.07 or -.3. As 7 
shifts discontinuously, a also makes a jump. In the direct scenario, a jumps from ^1.9 
to ~ 1.0 as 7 goes from positive to negative. For adjoint mixing, a goes from ^ 3.3 down 
to ~ 1.2. No longer seeing our (1/iV^) dependence, we need to explore the posteriors in 
more detail. 
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N 


I 


a 


direct 






3 


.716 


1.92 


4 


.602 


1.89 


5 


-.007 


.973 


oo 


-.069 


1.06 


V / / 1 oo 






adjoint 






2 


-.104 


.897 


3 


.706 


3.35 


4 


.594 


3.26 


5 


-.292 


1.22 


oo 


-.324 


1.29 


(l/^')|oo 






vertex 






3 


.512 


1.52 


4 


.472 


2.12 


5 


-.223 


1.28 


{l/N')\oo 







Table 7.2: Parameter values for (Teff for /3 = 0. 

Looking at the posterior for direct mixing first, Figures 17.31 through 17. 6[ we see that 
the peaks for N=3 and 4 are likely to be related, while starting at N=5 the peaks shifts 
dramatically. If we turn our attention to comparing the spectra, we can see what is 
driving the different peaks. The prime difficulty is getting the J=0 sector to fall into 
place. For N=4, Figure FTlj we see that the model is choosing to fit the conjugate splitting 
between the 0++ and the 0~~ at the expense of the magnitude of the radial excitations. 
(Such is also true for N=3.) When N=5, the model shifts to favoring matching the radial 
excitations at the expense of the low lying . Now the excited states are within ^ Icr, 
while the is about 2 ^ away. This behavior explains the shift in the peak of 
the posterior. What we would like is a model which can incorporate both the conjugate 
splitting as well as the radial excitations. The 3—2 states seem to be little affected by the 
shift of the peak - so far they are compatible with either emphasis. 

Turning to the indirect mechanisms, we see much the same behavior. For N—3 and 
4, we see one peak corresponding to fitting the conjugate splitting, while for N=5 and oo 
we see the other peak corresponding to fitting the radial excitations. Fitting the radial 
excitations gives a lower 7 than fitting the splitting, but the coefficient of elasticity is still 
positive, as the signal for zero elasticity is 7 ~ —1, from which the peak is still some 
distance away. 

Looking in detail at the spectrum resulting from the overall best fit (to the statistical 
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N = 3 
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Figure 7.3: 


The direct 


mixing posterior and spectrum for N=3. 



errors), the 0^^ notwithstanding, we see that the higher states are agreeing quite well 
with the hypothesis of spin misassignment. The M=4 states match the J=0 states with 
mixed quantum numbers quite well. When M=5, the ±+ state agrees with the J=l^+, 
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N = 4 



-1.5 




83 

O 
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O ^jx 



0.602 1.89 

++ -+ — +- 

Q I \ I \ I . L 

-1 1 2 3 4 5 

M = E m (n* + n ) 

Figure 7.4: The direct mixing posterior and spectrum for N=4. 



but the ±— state is predicted to be too heavy, yet the lattice data would agree will with 
the M==3. With the free predictions of the model giving such a signal, we should take 
seriously the issue of the spin assignments of the lattice operators. 
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N = 5 
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Figure 7.5: The direct mixing posterior and spectrum for N=5. 
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N 





M = Z m (n* + n , 

^ m m' 



Figure 7.6: The direct mixing posterior and spectrum for N=(X). 
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N = 3 
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Figure 7.7: The adjoint mixing posterior and spectrum for N=3. 
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N = 4 
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Figure 7.8: The adjoint mixing posterior and spectrum for N=4. 
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Figure 7.9: The adjoint mixing posterior and spectrum for N=5. 
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N 





Figure 7.10: The adjoint mixing posterior and spectrum for N=oo. 
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N = 3 
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Figure 7.11: The vertex mixing posterior and spectrum for N=3. 
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N = 4 
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Figure 7.12: The vertex mixing posterior and spectrum for N=4. 
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Figure 7.13: The vertex mixing posterior and spectrum for N=5. 
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N direct 


0++ 0~ 


2±+ 2±- 


3±+ 3±- 


4±+ 4±- 


1±+ i±- 


3 mjpc 
Tn*jpc 


4.308 6.228 
7.083 9.002 


6.838 8.758 
9.047 10.97 


8.005 9.924 
9.986 11.91 


9.108 11.03 
10.9 12.82 


10.15 12.07 
11.77 13.69 


4 


4.183 6.076 

6.994 8.887 


6.727 8.62 

8.965 10.86 


7.901 9.794 

9.907 11.8 


9.012 10.9 

10.82 12.71 


10.06 11.95 

11.7 13.6 


5 


4.276 5.249 
6.223 8.265 


6.9 7.873 
8.847 10.27 


8.114 9.087 
10.06 11.23 


9.263 10.24 
11.19 12.16 


10.35 11.32 
12.09 13.06 


oo 


4.105 5.164 
6.223 8.201 


6.737 7.796 
8.855 10.21 


7.955 9.014 
10.07 11.17 


9.108 10.17 
11.04 12.1 


10.2 11.26 
11.95 13.01 



Table 7.3: The spectrum for direct mixing with Ce// ^ at various N. For each N, the 
first mass is the lowest state in that channel, and the second mass is the radial excitation. 



N adjoint 


0++ 0~ 


2±+ 2±- 


3±+ 3±- 


4±+ 4±- 


1±+ 1±- 


3 rrijpc 
Tn*jpc 


4.25 6.215 
7.227 8.992 


6.963 8.746 
9.357 10.96 


8.214 9.913 
10.36 11.9 


9.396 11.02 
11.34 12.81 


10.51 12.06 
12.27 13.69 


4 


4.124 6.066 
7.134 8.879 


6.86 8.611 
9.281 10.85 


8.121 9.786 
10.29 11.79 


9.313 10.9 
11.27 12.71 


10.44 11.95 
12.21 13.59 


5 


4.178 4.855 
6.373 7.971 


7.021 7.516 

9.28 9.991 


8.317 8.748 
10.46 10.96 


9.534 9.916 
11.5 11.9 


10.68 11.02 

12.47 12.81 


oo 


4.06 4.81 
6.298 7.937 


6.925 7.475 
9.191 9.959 


8.229 8.709 
10.37 10.93 


9.455 9.879 
11.43 11.87 


10.61 10.98 
12.4 12.78 



Table 7.4: The spectrum for adjoint mixing with cTe// ^ at various N. For each N, the 
first mass is the lowest state in that channel, and the second mass is the radial excitation. 



N vertex 


0++ 0— 


2±+ 2±- 


3±+ 3±- 


4±+ 4±- 


1±+ i±- 


3 mjpc 
m jpc 


4.55 5.955 
6.973 8.796 


6.287 8.511 
8.565 10.77 


9.691 7.36 
11.72 9.5 


8.459 10.81 
10.43 12.63 


11.86 6.091 

13.52 8.793 


4 


4.33 5.901 
7.123 8.755 


6.301 8.463 
8.738 10.73 


9.645 7.439 
11.68 9.719 


8.592 10.76 
10.7 12.6 


11.82 6.612 
13.48 9.546 


4 


4.328 4.951 
6.782 8.042 


6.6 7.603 
8.52 10.06 


8.831 7.839 
11.03 9.725 


9.075 9.994 
10.91 11.96 


11.09 6.659 
12.87 10.22 



Table 7.5: The spectrum for vertex mixing with a^ff ^ at various N. For each N, the 
first mass is the lowest state in that channel, and the second mass is the radial excitation. 



Chapter 8 

Mixing With the K-string 



8.1 Adding another mixing parameter 

Recent investigations [43, 47] have lent support to the existence of a so-called k-string 
amongst the spectrum of pure gauge theory. The number and string tension of these new 
strings depends on N. For N=4 or 5, there exists a string with k=2. Following predictions 
from brane theory (see [47]), the string tension ak for a k-string goes as 



cTk _ sm ^ 



(8-1) 

CTi sm 

^k + 0{l/N^). (8.2) 



The k=:2 string may be schematically represented as the proximate union of two fundamen- 
tal strings, Figure 18.11 This bound state is lighter than two free flux tubes of equivalent 
length. As — > oo, the string tension approaches 2ai. For N=4, this ratio equals \/2, 
while for the physically interesting case of N=3, it equals unity. 

The important difference with the adjoint model is that this k-string carries an arrow 
of chromoelectric flux. Hence, eigenstates with negative charge conjugation may be built. 
This model possesses (at least) two towers of states: those arising from the traditional 
Isgur-Paton flux tube and states heavier in proportion by the ratio of the k-string tension 



(a) (Ti (b) o-j. 

Figure 8.1: The k-string. 
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to the fundamental sting tension. We are led to a basis consisting of four components: 
the L and R from the fundamental string, and another set of L and R from the k=2 
string. Eigenstates are built analogously between the two sectors, and so we expect the 
sub-Hamiltonians Hjp and Hk to be identical in form. With a four component basis, we 
write 

Hip M 
M Hk 



H = 



(8.3) 



where iJ/p and H^ are each two component direct mixing Hamiltonians with mixing ai 
(and elasticity 7). Hk is equivalent to Hip, Equation 14.21 with the replacement of a by 

Hk^Hip{a^ak). (8.4) 

The mixing matrix M is ours to determine (sensibly). Proposing that an L-orientated is 
much more likely to fluctuate into an L k-string rather than an R k-string (and vice versa), 
we write 

a2 
a2 



M 



(8.5) 



for some new mixing parameter 02- While we have now increased our parameter count by 
one, further analysis will yield quite an interesting result. 

Following earlier chapters, our merit function is defined and subject to minimization, 
only this time the parameter space is three-dimensional, making visualization of minimum 

surfaces slightly harder and the search for the global minimum more difficult. Note 
that two more states, the triply excited states at M — 0, have been added to the fits, 
making these not directly comparable to earlier ones. These states are in a range of 
mass compatible with J — A, but with congruent quantum numbers ++ and -, these 
states most likely have had their spins properly identified as 0. Pressing on, we map the 
posterior of a sensible range of parameter space and start the simplex from its best point, 
arriving at the spectra presented in Figures 18.21 through 18.81 and Tables 18.21 and 18.31 (An 
example of the posterior is shown in Figure [5.101 ) 

While the theoretical structure of the model carries for N down to 3 (or even 2), the 
calculations do not yield as pleasing a fit for such low N. At N=3, the ratio of the string 
tensions is unity, in which case we wonder if such a k-string really exists. The bases would 
be degenerate except for the identification of k-string states, and this is precisely how the 
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calculation is performed at = 3. However, the structure of the lattice spectra continues 
smoothly from high N down to N=2 [1], and so a model which performs well at higher N 
ought to remain applicable at lower N. Perhaps for N=3 there does exist a string heavier 
than the fundamental (such as the adjoint) but carrying an arrow of chromoelectric flux 
which no one has yet conjectured... Such speculations are beyond this work. 

For N=4 or 5, the results are much more pleasing, with the independence of the mixing 
parameters called into question by their incredible similarity, see Table I5TT1 Looking at the 
best fit parameters, we see that for these N the two mixing parameters are virtually equal! 
As the numerical calculation treats these numbers as independent, to see them converge 
might well be pointing at worthwhile physics to explore. At these N, this convergence 
of parameters lends some credibility to the model. We might even reconsider the k- 
string mixing model with the two mixing parameters replaced by a single parameter, thus 
reducing our parameter count again and allowing us to calculate some error bounds on 
our parameters. 



N 


7 


ai 


a2 




f3 = 










3 


1.49 


.925 


1.94 


81.8 


4 


.262 


.899 


.818 


34.5 


5 


.056 


.770 


.786 


27.9 


oo 


-.085 


.676 


.961 


14.3 


(3 = 0.05 










3 


1.45 


.848 


1.92 


7.44 


4 


.213 


.803 


.790 


6.80 


5 


.042 


.736 


.728 


7.36 


oo 


-.086 


.666 


.957 


5.59 



Table 8.1: Parameter values and for k-string mixing with 2 mixing parameters. 



85 



8.1 Adding another mixing parameter 



N k-string 


0++ 0— 


2±+ 2±- 


3±+ 3±- 


4±+ 4±- 


1±+ i±- 


3 m,jpc 
m jpc 


4.371 6.221 
6.908 8.758 


6.799 8.649 
8.829 10.68 


7.916 9.766 
9.745 11.60 


8.972 10.82 
10.63 12.48 


9.970 11.82 
11.49 13.34 


4 


4.290 6.088 

6.154 7.952 


6.950 8.749 

8.945 10.74 


8.175 9.974 

10.14 11.94 


9.332 11.13 

11.17 12.97 


10.42 12.22 

12.11 13.91 


5 


4.218 5.757 
6.098 7.637 


6.911 8.450 
9.010 10.55 


8.149 9.688 
10.19 11.73 


9.318 10.86 
11.20 12.74 


10.42 11.96 
12.14 13.68 


oo 


4.081 5.432 
6.279 7.853 


6.817 8.168 
9.155 10.29 


8.071 9.422 
10.26 11.39 


9.252 10.60 
11.28 12.42 


10.37 11.72 
12.24 13.38 



Table 8.2: The spectrum for k-string mixing at various N. For each N, the first mass is 
the lowest state in that channel, and the second mass is the radial excitation. 



N k-string 


0++ 0~ 


2±+ 2±- 


3±+ 3±- 


4±+ 4±- 


1±+ i±- 


3 rrijpc 
Tn*jpc 


4.413 6.109 
6.963 8.659 


6.847 8.542 
8.885 10.58 


7.967 9.662 
9.803 11.50 


9.025 10.72 
10.69 12.39 


10.03 11.72 
11.55 13.24 


4 


4.334 5.939 
6.134 7.740 


7.001 8.606 

8.950 10.56 


8.229 9.834 
10.17 11.78 


9.389 10.99 
11.22 12.83 


10.48 12.09 
12.17 13.78 


5 


4.271 5.743 
6.068 7.540 


6.961 8.433 
9.020 10.49 


8.197 9.669 
10.23 11.70 


9.364 10.84 
11.25 12.72 


10.47 11.94 
12.19 13.66 


oo 


4.093 5.425 
6.242 7.575 


6.828 8.161 
9.112 10.45 


8.082 9.414 
10.21 11.54 


9.263 10.59 
11.20 12.53 


10.38 11.71 
12.13 13.46 



Tabic 8.3: The spectrum for k-string mixing with adjusted weightings at various N. For 
each N, the first mass is the lowest state in that channel, and the second mass is the radial 
excitation. 



86 



8.1 Adding another mixing parameter 



95.2 



81.6 



N = 3 



t X I 



X O 



I 



-rO 



1 



I 



X 



O 



1 .49 0.925 1 .94 



++ -+ — +- 



b0a2newkN3.mat 



M = E m (n"^ + n ) 

^ m m' 



Figure 8.2: The spectrum for N=3. 
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Figure 8.3: The spectrum for N=4. 
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52 



27.9 



N = 5 
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Figure 8.4: The spectrum for N=5. 
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32.1 



14.3 



I 



++ -+ — +- 



tO 



I 



-0.0854 0.676 0.961 - 



b0a2newkNoo.mat 
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Figure 8.5: The spectrum for N=oo. 
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9.58 



7.44 



^ o 



N = 3 



O 



1 .45 0.848 1 .92 



++ -+ — +- 



b5a2newkN3.mat 
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^ m m' 



Figure 8.6: The spectrum for N=3 with adjusted weightings. 
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15.5 
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Figure 8.7: The spectrum for N=4 with adjusted weightings. 
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12.8 



7.36 



I 



N = 5 



0.0417 0.736 0.728 
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Figure 8.8: The spectrum for N=5 with adjusted weightings. 
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9.56 



5.59 
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-0.0858 0.666 0.957 - 
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Figure 8.9: The spectrum for N=oo with adjusted weightings. 
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Figure 8.10: The posterior for N=5 with adjusted weightings 
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8.2 Combining the mixing parameters 

Returning our attention to Equation 18.21 we now suppose that the mixing matrix M 
does not contain a new parameter 02 but instead uses our old parameter ai to perform the 
mixing between sectors. That such an identification of parameters may be performed is not 
immediately obvious: ai stems from the direct mixing between the L and R orientations 
of the same type of string, while 02 represents mixing between the two types of string. 
However, the results immediately above strongly suggest that, at least for > 4, such 
identification may be valid. As always, the proof is in the pudding. 

We repeat the calculations (again), but this time only for N — 4, 5,andcx). In Ta- 
ble [5]3] we present the resulting best fit parameters and their associated x^- The details 
of the spectra are given in Figures 18.111 through 18.151 and Tables 18.51 and 18.61 The most 
interesting feature of this model is that it works, and quite well at that. The values of 

are comparable or even better than the earlier models. However, we must not put 
too much faith in the k-string model as a successful explanation of the lattice spectrum, 
as the continuation down to N=2 which manifests in the lattice data cannot be repro- 
duced. Nonetheless, the k-string model must have some validity: it produces (for adjusted 
weightings) a x V d.o.f. of 7 / (8 - 2) « 1! 

Ultimately, the complete spectrum of pure gauge theory is unlikely to be adequately 
described by such simple effective Hamiltonians as these. However, by distilling the mul- 
titudinous gauge degrees of freedom down to a few clearly defined parameters, we hope to 
capture the essential aspects necessary to describe the gluonic physics relevant to the real 
world. 



N 


7 


(cr^) 


a 






/3 = 












4 


.276 


(6) 


.878 


(26) 


35.4 


5 


.054 


(38) 


.774 


(39) 


27.9 


00 


-.080 


(45) 


.779 


(49) 


17.2 


(3 = .05 












4 


.215 


(1) 


.799 


(1) 


6.81 


5 


.043 


(76) 


.735 


(92) 


7.36 


00 


-.112 


(79) 


.720 


(100) 


6.31 



Table 8.4: Parameter values and for k-string mixing with 1 mixing parameter. 
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Figure 8.11: The spectrum for N=4. 
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Figure 8.12: The spectrum for N=5. 
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16r 
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Figure 8.13: The spectrum for N=oo. 
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Figure 8.14: The spectrum for N=4 with adjusted weightings. 
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Figure 8.15: The spectrum for N=5 with adjusted weightings. 
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Figure 8.16: The spectrum for N=oo with adjusted weightings. 
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state 


SU(4) 


SU(5) 


SU(oo) 


0++ 


4.283 


4.219 


4.105 


0++* 


6.24 


6.08 


6.10 


Q+ + ** 


7.40 


7.43 


7.39 


2±+ 


6.95 


6.91 


6.81 


2±+* 


8.99 


9.00 


9.09 


4-+ 


9.33 


9.32 


9.23 


4-+* 


11.18 


11.20 


11.16 


1±+ 


10.42 


10.42 


10.34 


]^±+* 


12.11 


12.14 


12.01 


0~ 


6.04 


5.77 


5.66 


Q — * 


8.00 


7.63 


7.66 


Q 


9.15 


8.98 


8.95 


2±- 


8.70 


8.46 


8.37 


2±-* 


10.75 


10.55 


10.65 


4+- 


11.09 


10.87 


10.79 


4+-* 


12.93 


12.75 


12.72 


1±- 


12.18 


11.97 


11.89 




13.87 


13.69 


13.65 



Table 8.5: The spectrum for k-string mixing with 1 parameter. 



St at (> 


SIT(4) 


SU(5) 




0++ 


4.332 


4.270 


4.148 


0++* 


6.15 


6.07 


6.04 


Q + + ** 


7.45 


7.46 


7.41 


2±+ 


7.00 


6.96 


6.86 


2±+* 


8.96 


9.02 


9.11 


4-+ 


9.39 


9.36 


9.27 


4-+* 


11.23 


11.25 


11.21 


1±+ 


10.48 


10.4G 


10.38 


1 1 


12.17 


12.12 


12.13 


0— 


5.93 


5.74 


5.59 


Q — * 


7.74 


7.54 


7.48 


Q 


9.05 


8.93 


8.85 


2±- 


8.60 


8.43 


8.30 


2±-* 


10.56 


10.49 


10.55 


4+- 


10.99 


10.83 


10.71 


4+-* 


12.82 


12.72 


12.64 


1±- 


12.08 


11.93 


11.82 


1 1 


13.77 


13.(i(i 


13.57 



Table 8.6: The spectrum for k-string mixing with 1 parameter with adjusted weightings. 



Chapter 9 

A Novel Construction of 
Lattice Operators 

Working with the Isgur-Paton model has brought to light the importance of correctly 
identifying the spin of states measured by lattice operators. For a review of lattice gauge 
theory, sec [4,45,48]. While the lattice is an excellent method for determining the gross 
properties of a state, such as its mass or branching ratios, it can give only limited infor- 
mation regarding a state's structure. Yet that structure has a strong impact on the state's 
quantum numbers. The modern technology of smeared lattice operators [49] further com- 
plicates extracting structure information from correlation functions. As we have already 
seen, differentiating between models of these gluonic states requires understanding the 
spin structure of the entire spectrum. 

Our problem demonstrates that what is a solution to one problem can become another 
problem in itself. In the continuum, field theory is plagued by ultraviolet divergences [50]. 
Introducing the lattice discretization of space-time successfully regularizes these diver- 
gences, but at a cost. By discretizing space-time, we break the continuous rotational 
symmetry down to that of the lattice, so that rotations of arbitrary degree are reduced to 
a discrete set of rotations. For conventional square lattices, the rotation group is broken 
to the discrete group 0(2), rotations by multiples of 7r/2. 

Let R{9) be the operator for a rotation of a loop by an angle 6. Then, its action on a 
loop may be written 

i;e = R{e)i^e=o- (9.1) 
in.-? 
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To represent a state of arbitrary J in the continuum, we then write 

*j = / d^e^'^Ve) |0 >, (9.2) 
Jo 

where ipe^o represents the state along some fixed axes, which is then multiphed by the 
phase iJ9 when rotated to angle 9, and |0 > is the vacuum. On the lattice, this becomes 

3 3 

VI/,, ^ ^ e'''-Mo.) |0 >- |0 > . (9.3) 

n=0 n=0 

Accordingly, the lattice can only distinguish rotations of tt/2. This rotational ambiguity 
shows up as an ambiguity in the spin represented by the state. For example, taking J=4, 

3 

*j=4 = ^e™2>fl(e„) |0> (9.4) 

3 

= E^«(«")|0> (9-5) 
= *j=o, (9.6) 

thus an operator which couples to J=0 will also couple to 3—4 on a cubic lattice. So far, 
lattice methods only extract the two or three lowest states in a given symmetry channel, 
and the convention of assigning the lowest compatible spin seems applicable, but as we 
continue to explore the heavier gluonic states, we need to critically reevaluate our methods. 

Earlier work [44] explored a technique to distinguishing the J=4 and the J=0 in the 
PC = ++ sector on D=2+l lattices, but its result was inconclusive. Recently progress 
has been made in distinguishing the J=4 from the J=0 on D=3+l cubic lattices [51]. 
By methodically and painstakenly sorting out the appropriate symmetry contributions for 
quite twisted (and smeared) conventional lattice operators, those authors claim to separate 
out the J=4 channel from the J=0. However, their technique is quite specific to the spins 
and 4 and the 3+1 cubic lattice. What we are after is a technique that may be applied 
to arbitrary spins in arbitrary dimensions (and indeed even on arbitrary lattices). 

To overcome our difficulties with arbitrary spin resolution on a cubic lattice, we need 
operators which we can place onto the lattice at arbitrary angles with a minimum of 
distortion. These operators may be placed on the lattice at relative angles other than 
7r/2, thereby picking up different phase contributions than operators constrained to align 

^For the initial idea to try this new method, one must thank Mike Teper. 
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to the axes of the lattice. These phase contributions will not be identical for spins which 
are equal moduloA, and thus the J=4 will be distinguishable from the J=:0. Such new 
operators must still be smooth on physical length scales, requiring the development of the 
arbitrary smeared link between any two spatial lattice sites. 

Conventional operators are constructed from closed contours of gluonic links. The 
requirement of gauge invariance is what necessitates using closed contours; otherwise under 
a gauge transformation U{x) V {x)U fi{x)V^ {x + fi) the operator acquires unmatched 
contributions from its endpoints which break gauge invariance. Thus, the simplest gluonic 
operator is the real part of the trace of the basic plaquette, which has the quantum 
numbers of the vacuum. To find the C = — analogue of the operator we take the imaginary 
part of the trace. As the plaquette is rotationally invariant, to construct an operator of 
nonzero spin, we must use operators with more interesting shapes. States with i—2 can 



be made from operators which are longer than they are broad, Figure 9.1(a) To couple 



to J=l, operators with very little intrinsic symmetry must be used, eg Figure 9.1(b) To 
construct an operator of negative parity, basic shapes of negative parity are combined in 
an appropriate linear combination. As we go to the continuum limit, a — > 0, the size of our 
basic operator becomes infinitesimal, while the size of the physical giueball remains finite. 
To overcome this difficulty, basic links are replaced by their smeared counterparts [49]. 
Smearing increases both the length and breadth of the links used. As the level of iterative 
smearing increases, the resultant smeared link is a complicated combination of the raw 
links of the lattice which is "smooth" on physical length scales. These techniques, while 
improving the overlap with physical states, prevent an easy visualization of the structure 
of the state. 

Adopting a broader perspective, we would like to construct operators using an arbitrary 



(a) Operators coupling to J=2 (b) Operators coupling to J=l 

Figure 9.1: Conventional lattice operators. 
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selection of lattice points. What we need is a "smeared link" between any two points in a 
timeslice, not just along an axis. Labelling the link starting at i between adjacent lattice 
sites i and j by Uij , we construct the matrix M for each timeslice by 



Uij if Uij exists 
otherwise 



(9.7) 



M is an n X n sparse matrix with a high bandwidth. As Uji — C//^ , Mji — Mj-, and M is 
Hermitian. Using the well known expansion, we define 

K = (9.8) 

= 1 + aM + a^M^ + . . . (9.9) 

with parameter a. When M is multiplied by itself, nonzero entries are the sum of the 
products of the link variables between the sites consisting of a number of links equal to 
the number of times M appears. The full matrix K represents a matrix of propagators 
between every site in our timeslice, i.e. Kij = Kj^ is the propagator from site i to site 
j. Within its radius of convergence, the expansion for K can be used to determine which 
paths are contributing to K. The weighting parameter a controls the extent to which 
longer paths are contributing. For very small a, 

K Kil + aM + aM'^. (9.10) 

Looking at the main diagonal of K, paths of length 2 are contributing, thus 

if,, wl + 4a2, (9.11) 

as there are 4 paths from the center site i to its nearest neighbors and back again. Fig- 
ure 19.21 These paths contribute unity each as the product along a link and back again is 
simply U X W. For slightly greater a, along the diagonal 

Ku^l + o?[M'% + a^[M%,. (9.12) 



Now that we have paths of length 4, there are 12 paths which contribute unity each 
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Figure 9.2: Four paths of length 2. 



to [Af^Jii, three off of each link coming from i, Figure 9.3(a) And we also have the 4 



plaquettes in the plane which have site i as a vertex, Figure 9.3(b) To illustrate our 



(a) Twelve paths of (b) Four plaquettes 

unity 

Figure 9.3: Paths of length 4. 
method, we can estimate the expectation value of the plaquette by 

< \] >Ki — (Tr K - 1 - Aa^ - 12a^) (9.13) 
An 

As a increases, paths of increasing length contribute to K. We must beware, however, 
of the divergence of K as a approaches some critical value ac which makes the denominator 
of Equation (|9.9p equal zero. To find this critical value, we can expand the trace of K as 
a sum over terms containing the reciprocal of the eigenvalues of M: 

TrK = Tr{l + aM + a^M'^ + ...) (9.14) 
= Trl + aTrM + a'^TrM^ + ... (9.15) 
= n + a^Ai +a2^A2 + ... (9.16) 

n 

= ^{l + aX,+a^X^^ + ...) (9.17) 

i=l 
n 

2—1 

which can equal oo if any Ai equals a. As we take a as ranging from to ac, it is the 
largest eigenvalue of M which determines the divergence of Tr K . A plot of < Tr K > 
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is shown in Figure 111.11 - we can easily see the divergence at ac as weh as the change in 
sign when a > ac- To be effective in constructing operators, we need a large enough that 
paths significantly longer than the shortest path between two sites i and j contribute to 
Kij, while not so large that we are diverging. The range of effective a in this calculation 
is found empirically. 

From our matrix K, we can construct an operator by taking a closed path of propaga- 
tors, so that, e.g. 

O = K,,K,kKk^ (9.19) 

defines a triangular operator. Figure 19.41 While triangular operators are quite pretty, 
they possess too much symmetry, coupling effectively only to J=0. In order to have less 



k 




Figure 9.4: Triangle based lattice operators. 

symmetry, we construct operators by joining two triangles together to produce a rhombus, 
Figure 19.51 Each rhombus can be put onto the lattice in one of three orientations, but 
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Figure 9.5: Rhombus based lattice operators. 

as we are still working with a cubic lattice, these orientations will have slightly different 
actual sizes, and hence different vacuum expectation values. (Actually, the two "vertical" 
orientations should be equivalent, while the "horizontal" orientation will be different.) A 
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solution, following [44] , is to normalize the operators so that 

<Ol >^cl^<Ol > forj^l. (9.20) 

Such normalization should account for any directional variance. From these normalized 
operators we can construct our physical operators by multiplying each base operator by 
a phase iJ9 corresponding to its desired spin and its angle of orientation. We take the 
three orientations shown in Figure [9751 as well as their 7r/2 counterparts as our physical 
operators, for a total of six base operators. We need pairs of loops with a relative angle 
of 7r/2 to produce good 3—2 operators as well as to affect the necessary cancellations 
detailed below. Next we take the zero momentum sum of O across each timeslice. From 
these numbers we can calculate the correlation function 

Co{t) =< O{t)\O{0) > (9.21) 

which decomposes into a sum of decaying exponentials 

Co(t) = ^a„e-™"*. (9.22) 

n 

For large t, we can extract the lowest mass from 

Coit) oc e-™°(*) + e-"<'(^-*) (9.23) 

in Euclidean space-time, where the second term arises from our finite lattice with time 
extent T=L. As we are interested in distinguishing the 4++ from the 0++, as well as from 
the vacuum which has the same quantum numbers, in practice we use vacuum subtracted 
correlation functions 

Co{t) Co{t)- < O >^ . (9.24) 

Operator construction is as much an art as it is science. Finding the right combination 
of shapes and smearing level which couple to the desired state is educated guesswork. 
Our first operators are based on two-triangle rhombi, which may be rotated by 7r/3 while 
remaining roughly the same size. Figure 19.51 As these operators are at relative angles 
other than 7r/2, the J=0 and J=4 will now have different phases on each rhombus, and 
hence are distinguishable. 
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(a) Unwanted torelon contribu- 
tion. 



-liV 



(b) Removing torelon contribu- 
tions. 



Figure 9.6: Torelon contributions on the lattice. 



One problem encountered is when flux from torelons contributes to Kij. This flux 
becomes a problem at the higher a necessary to get a significant contribution to K. As a 
increases, not only flux along the shortest path between sites i and j but also flux stretching 



the long way around the lattice will contribute, Figure 9.6(a) As we are interested in 
operators based on flux along the shortest closed path, this extra flux needs to be removed. 
Luckily, a technique for removing torelon flux contributions exists. For each timeslice, 
the links perpendicular to a boundary layer in either direction are multiplied by -1 (see 



Figure 9.6(b) ) and the calculation of the operators is repeated. When the values are 
combined with the earlier values, flux which crosses a boundary only once will cancel 
away, leaving the flux contributions that we desire. In practice, this procedure is done 
three times: once each for either direction alone, and once for both directions together. 
Thus we increase our computational load by a factor of 4, not an insignificant amount. 

Operators of various sizes may be constructed by these means. We label the size of 
the operator by half the number of lattice units on a side of its constituent equilateral 
triangles, so that, eg. Figure is a "size 2" triangle, and the full operator Figure [^31 
is also a "size 2" . We should note that these are not perfect equilateral triangles, owing 
to the minor distortion necessary to put the vertices onto sites of the lattice, and so the 
angles involved are not perfect multiples of 7r/3. A simple center-of-mass calculation for 
the loops yields the true angles, but in practice the difference is miniscule, and for these 
calculations perfect 7r/3 rotations are assumed. 

Let's look in detail at why we need to use two sets of rhombi rotated by tt/2. Suppose 
we have loops Ij for j — I...N identical (in the continuum) upto rotations of 6j, where 
J2j (^j = 27r(mod(?i). (The mod^i simply accounts for whatever offset the first loop is at.) 
For example, the operator in Figure [9751 has 6j E {7r/6, 57r/6, 37r/2}. Consider as a trial 
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operator for definite spin J 

Lj = ^e*^^^7j. (9.25) 

j 

Does this mean that its projection onto the vacuum reaUy is of spin J, ie 

{n\Lj\J) (xSjj7 (9.26) 

Consider 

(4L,) = 5^e^(«^-«.')^(/,./,). (9.27) 

jj' 

Inserting a complete set of states J^j i^^o the expectation value on the r.h.s, 

{i,d,) = Y,mj'\j){j\m (9.28) 

by rotating Iji over to Ij, 

J 

(9.30) 

J 

Thus, 

(4L,)^^Cj^e^(«.'-^^)(^-^~) (9.31) 
J Jj' 

= Y,CjF{J,J), (9.32) 
J 

where 

FiJJ) ^ (9.33) 
j j' 

Thus our state with spin J will have no overlap with states of spin J as long as either 
factor of F{JJ) is zero. 

For example, let's consider using only the three loops shown in Figure 19.51 as our 
operator, with J=4 and 9j £ {7r/6, 57r/6, 37r/2}. Evaluating the first factor of F{JJ) we 
have: 

J = : e-^'i + e**T + e^'^ = (cool...) (9.34) 
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J = 2 : e^^f + e^'T + g^'* = o (cool...) (9.35) 

J = 3 : e^'t + e^^if + e^*^ = (cool...) (9.36) 
J = 1 : e^^t + e^** + e^'^ = 3i (ouch!) (9.37) 

And so this operator is going to have some overlap with the J=l state. Now consider 
adding three more loops to the operator with angles 9j G {7r/2, 77r/6, ll7r/6}. Looking 
only at its overlap with J = 1, 

J = 1 : + e'^^ + e'^^ = -3i (9.38) 

which cancels the earlier contribution of 3i from J = 1. Thus, by using six loops at the 
angles given above, we can construct an operator with J=4 and no overlap with states 
having J < 3. 



Chapter 10 



Performing the Calculation 



For this calculation wc work with SU(2) on a D=2+l Euclidean lattice with L=16 and 
/? = 6. As we are in pure gauge theory, the only field present is the gauge field defined on 
the links between lattice sites: 



where Xi gives the site where the link begins and /i is a unit vector in the direction of the 
link. Working with the gauge group SU(2), the link variables are SU(2) matrices, which 
require four real numbers to represent. The gauge field is updated using the Kennedy- 
Pendleton heatbath algorithm [52] supplemented by over-relaxation sweeps [53] and the 
occasional global gauge transformation. Calculations are performed on a stored set of 
1000 configurations separated by 40 heatbath sweeps each. Since we want to calculate the 
rest energy, we take the zero momentum average of our operators in each timeslice. 



One Euclidean direction is taken to be the time direction, say, and sites with a common 
value of comprise a single timeslice. There are a total of lattice sites, and it is 
convenient to define n = to be the number of sites in a given timeslice. 



U^{xi) = exp / A{x)dx, 



(10.1) 




(10.2) 
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10.1 Building M 

First, a set of configurations is generated using standard techniques. After an initial 
period of 1000 sweeps thermalizing the lattice, configurations are saved every 40 sweeps. 
(Actually, they were saved after every 10 sweeps, but to reduce the computational time to 
something reasonable, only one quarter of that data was used.) Measurements are then 
performed on each configuration as follows. 

For each timeslice, sites are labelled by an index i — l...n. Then, the n x n matrix M 
is formed from the gauge fields U such that the entry for the link from site i to site j is 

M,, = = C/^,. (10.3) 

Note that AI = A/^ = 1, so that M is Hermitian. Because the gauge fields U have a 
nearest-neighbor relationship, only certain entries of M are non-zero. For each row i of M, 
the only non-zero entries are in the four columns j corresponding to the four neighbors of 
i. This construction leads to a very sparse matrix M having a band diagonal structure. 
For L=16, the diagonals of M are located at offsets (where means the main diagonal): 

diags = ±1, 2, 3, 29, 30, 31, 32, 33, 479, 480, 481, (10.4) 

See Figure 
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Figure 10.1: The band structure of M. 
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10.2 Hermitian Inverse by LU Decomposition 



We wish to calculate 



and so M is multiplied by —a and added to the n x n identity matrix to form K~^. 
Next comes the most computationally expensive part of the calculation: inverting K~^. 
Several algorithms were tested for speed and accuracy: the one used in the calculation 
was a modified version of the Super LU program, version 1.1 [54]. SuperLU is a good tool 
to use when inverting any general sparse matrix, however an algorithm specific to the 
structure at hand will prove useful in future calculations. 



10.2 Hermitian Inverse by LU Decomposition 

While the following algorithm was not implemented in the current work because of time 
restrictions, any future work on this technique will find it usefuU. The high bandwidth 
of means there is much filling in of non-zero entries as the inverse is calculated. 
The following algorithm seeks to minimize such fill-ins as well as storage requirements, 
necessary steps when building a production version of this calculation to handle lattices 
larger than L=16. First, we write 



K-^ = l-aM 

= LxU, 



(10.6) 

(10.7) 



where K ^ is Hermitian and L and U are complex. Taking the 3x3 case as a graphical 
example, let's display all the degrees-of-freedom: 
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Counting degrees of freedom: 



-real - 



imag n{n + 1) ® n(n + 1) 



(lO.i 



(10.9) 



The l.h.s has degrees of freedom, while the r.h.s. has 2 x n{n + 1). If we impose the 



conditions U = and La = L*^, then we are back down to on the r.h.s. As the entries 



116 



10.2 Hermitian Inverse by LU Decomposition 



on the diagonal are unity while all other entries have a magnitude less than unity, the 
property of diagonal dominance insures we can ignore pivoting. The storage requirements 
for the algorithm can be met with one full, complex matrix of size n x n. First, store L 
in the lower part of that matrix via forward substitution (read Numerical Recipes [34] for 
more details of forward and backward substitution): 



V i<j 



1 

Lij = -r^{K-/-Y,LikL%) fovi>j (10.10) 



fe=i 



Then compute L~^, using a transposed storage scheme into the upper portion of our 
storage matrix and neglecting to store the diagonal elements (which are trivially related 
to L's diagonal elements): 



^7/ = 7^ ; L-i' = -F^y.L.kL^^ forz>j 



(10.11) 



Next, compute K, where the tilde indicates that we are close but not quite there yet, using 
back substitution and storing over L~^: 

T-l 1 " 

If we write K = K + 6, some iterative improvement gives us the error 5 



k 

-1 

Sij = -Z ^ ^ -^ik^kj 
^" fe=l 



Sn = J 



(10.13) 
(10.14) 



33 



which may be subtracted off from K to give our final K. Schematically, the storage re- 
quirements go like: 











■•. K 




'■■ k 




■•. k 


L ■•• 




L 




L '■■ 




A '■■ 




6 ■■. 



until we fill the entire matrix with K. The last stage of iterative improvement may be done 
efficiently by taking into account the sparsity of K~^. 
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10.3 Evaluating the operators 



10.3 Evaluating the operators 

Once we have our matrix K we can begin evaluating our operators. We use our rhomboid 
operators at sizes 3,4,5, and 6. Size 2 was also computed, but never showed much promise 
for being useful, and so is hereafter neglected. For each site in a timeslice we compute 
the values for all six loops and all four sizes (and for several values of alpha). Then we 
compute the normalizing constants ce for the loops via Equation 19. 201 If everything were 
perfectly symmetric, these would all be unity, but the slight distortions caused by the 
cubic lattice will produce non-zero normalizing constants. A measure of how similar our 
loops are is given by cg. We now have the choice of applying our constants cq either 
to all of the loops at hand or only to those two loops which are obviously going to be 
different than the other four, namely loop 3 in Figure WIE\ and its tt/2 analogue. If we 
apply our constants to the two loops only, we are expecting the tt/2 analogues to produce 
cancellations (over the duration of the calculation) when computing the J— 2 so that there 
is no overlap with the J=0. If, however, we apply our constants to all the loops, we expect 
no fortuitous cancellations and instead rely on the actual phases applied to our loops to 
affect the cancellations. In practice, there is not much difference between the two methods. 
If anything, the more general method of normalizing all the loops does produce a slightly 
better result and will be the focus for the remainder of this work. We can measure this 
effect by looking at the overlap of our J=2 operator with the J=0: ideally, subtracting 7r/2 
analogues directly, ie 4'J=2 = ''Po — ^e+Ti/2i will produce an overlap equal to the statistical 
error of the calculation. As we are looking for a method which can be applied to even 
more general situations than the current one, we normalize all the loops and believe the 
calculation of the overlap to be the true overlap of our operators (upto statistical noise). 

We now have a set of operators O with several indices: size, loop orientation, and a. 
Actually, we have four sets of these operators in order to cancel the torelons. An important 
detail is that the torelon cancellation must be done before the normalizing constants are 
applied. In order for the torelon flux to cancel, the four sets of loops are averaged together 
to produce a torelon canceled set of loops, and then the normalizing constants applied. 
Any other ordering of these steps will not work. We wish to reduce all these numbers down 
to a few correlation functions from which we can estimate masses. To produce an operator 
coupling to a particular spin J, we add our six loops together with the appropriate phases. 
Note that the center of the first loop, number 1 in Figure is at angle 7r/6. Including 
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10.3 Evaluating the operators 



notation for the normalizing constants, we have explicitly: 

Oj=o = tp7r/6 + C5,r/6V'57r/6 + C3,r/2V'37r/2 + 

Cn/2i'7r/2 + C77r/6V'77r/6 + Cii7r/6 V'Htt/B (10.15) 

e'*"/'c,/2V'V2 + e'''"/'c7V6V'7./6 + e''""/'cnv6V'iiV6 (10-16) 

e*-/2c,/2V'./2 + e^'^'^/<5c7V6V'7V6 + e''""/'cn./6V'ii,r/6 (10.17) 

As the vacuum has the same quantum numbers as our ground state, 0++, we need to use 
vacuum-subtracted operators for the J=0 correlation function, and as we are unsure of 
what is happening at J=4, we use vacuum-subtracted operators there too, just for good 
measure. 

= 0- <0> (10.18) 

Now that we have our operators calculated for each timeslice, we are in position to 
compute our correlation functions: 

Cj{t) = < 6/it)\dj{o) > (10.19) 

= < {o*j(ty < o*j >)\{Oj(oy <Oj>)> (10.20) 

= < O''j{t)\Oj{0) > - <0} xOj > (10.21) 

These correlation function will be normalized such that Cj(0) = 1, ie 

Cj{t) ^ Cj{t)/Cj{0) (10.22) 

Similarly, we can define and compute the normalized cross-correlation for states of different 
J, eg: 

C^t) = <Om\Oo> (10.23) 



^J< 0^(^)104(0) X dm\Oo{o) 



> 
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10.4 Error analysis 



If we take t = in Equation ll0.23[ we have computed the actual overlap of our operators. 
This number will be our measure of how well we are distinguishing the J=:4 from the J=0. 
From these correlation functions we can extract effective masses in units of the lattice 
spacing a: 

</W = ln^^ (10.24) 

10.4 Error analysis 

To estimate the errors on our secondary quantities, ie our correlation functions C and 
effective masses m, we perform a standard jackknife analysis [48]. Consider a primary 
quantity x with N measured values Xi,i — 1...N. The sample average is simply x = 
Now consider a secondary quantity y = y{x). Our best estimate of y is given 

by 



y = y{x) + y{x) (10.25) 

The jackknife procedure estimates the error on y by evaluating the distribution of N 
jackknife estimates y'^ calculated with one sample missing from x for each, ie: 



^i^T^y.^i (10.26) 



S TV - 1 ^ ■ 



Then, the jackknife estimates are yj' — y{Xj) with mean 

1 ^ 



N 



and variance 

4' = ^hyi~y')' (10.28) 

For primary quantities this variance is equivalent to the simple variance of the measure- 
ments, while for secondary quantities (such as effective masses) it produces a more stable 
estimate of the error 

y = y±ay., (10.29) 
These are the errors reported in the following chapter. 



Chapter 11 

Results From the Lattice 



11.1 Preliminaries 



First, let's look at the behavior of < Trif > as a function of a. As < TrK >, ioi 
small enough a, is a measure of the expectation value of a plaquette, upto some constants, 
its value gives us an idea of how our calculation is performing. As our matrix M is in 
the denominator of our expression for K, Equation 19.91 there exists a critical ac such 
that K diverges. We can estimate this ac by evaluating the Frobenius norm of M. In 
general, M will be a complex matrix, and its norm will be difficult to evaluate analytically. 
However, we can derive our result by considering a simple limiting case: that of a minimal 
two dimensional lattice at its "cold" start before thermalization. This minimal lattice 
possesses 4 sites arranged as a square (or plaquette). Using the symbols 





1 







V = 




and C = 






1 








(11.1) 



M is given by 

C '7 '7 C 

V C C V 

V C C V 
C V V C 

The Frobenius norm of this matrix is easily calculated to equal 4. Thus, with an ac — l/A 
the norm of the denominator of Equation l9.9l will equal 0, and K^^ will be non-invertible. 



Mc, 



(11.2) 
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11.1 Preliminaries 



Tr(K) vsa 

300 I 1 1 




0.05 0.1 0.15 0.2 0.25 0.3 

a 

Figure 11.1: < Tr K > as a. function of a. 

For a "hot" configuration, the unitary matrices making up M wiU be rotated away from 
their principle axes, and the norm of M wiU be shghtly less than 4, so that ac > -25. 
While these are heuristic arguments, we see precisely this behavior in actual calculation. 
As a 0, < Tr K > approaches < [] >; as a ^ .25, < Tr K > diverges. Figure [11.11 
plots Eauation l9.13l as a function of a. We easily see < Tr K > diverging as it approaches 
some a ~ .25. Actually, the final positive point is at a = .26, indicating that ac > .25. 
The last point, at a — .3, has a negative value for < Tr K >, an unphysical result if 
< Tr K > is to represent the expectation value of some sum of closed loops on the lattice. 
Calculations performed at a > ac yield operators with a divergent vacuum expectation 
value and correlation functions which fall to after only one lattice spacing. While this 
regime is not useful to the present calculation, perhaps future work will find something 
worthwhile to explore. 

Next, let's look at the normalizing constants ce which are applied to our individual 
loops to overcome any spatial rotational discrepancies our finite cubic lattice imposes. In 
Table [TTTT] we display these constants for the six loops at the four values of a used in this 
calculation. We notice that as a increases these constants approach unity. Apparently, 
as longer and more spread-out paths contribute to the loop, the orientation of the loop 
on the lattice becomes less significant, which makes sense. The size 6 loops show the 
most variation between orientations, to the order of 10'' at a = .1! Even at the more 
modest a — .2 the variation is of the order lO'^. If we did not apply these normalization 
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11.2 J = and the effect of torelons 



constants, we would not get a balanced contribution from the different orientations of loop, 
but perhaps we should be concerned about using such large normalization constants! At 
least the behaviour of these constants is doing what we would expect: the four equivalent 
orientations have constants consistently close to unity, the smallest loop needs the least 
correction, and the largest loop has its smallest (non-unity) constants close to ac where 
the extreme spread of the loops should help cancel the difference in orientation. Thus, we 
proceed to take these constants at face value and implement them in the calculation. 



loop 




size 




# 


3 


4 


5 


6 






a = 


.1 




1 


l.OOOOe+00 


l.OOOOe+00 


l.OOOOe+00 


l.OOOOe+00 


2 


1.0082e+00 


9.6166e-01 


1.0237e+00 


1.0161e+00 


3 


1.8708e+01 


1.3566e+02 


1.3982e+05 


1.4130e+07 


4 


1.0055e+00 


9.5868e-01 


1.0346e+00 


1.0211e+00 


5 


9.9742e-01 


9.7717e-01 


1.0386e+00 


1.0196e+00 


6 


1.8997e+01 


1.4096e+02 


1.4026e+05 


1.4339e+07 






a = 


.15 




1 


l.OOOOe+00 


l.OOOOe+00 


l.OOOOe+00 


l.OOOOe+00 


2 


9.9819e-01 


9.5151e-01 


1.0282e+00 


1.0201e+00 


3 


4.6171e+00 


1.7857e+01 


2.5579e+03 


8.5268e+04 


4 


1.0044e+00 


9.5565e-01 


1.0428e+00 


1.0210e+00 


5 


9.9262e-01 


9.7976e-01 


1.0322e+00 


1.0253e+00 


6 


4.7088e+00 


1.8702e+01 


2.5747e+03 


8.6693e+04 






a. = 


.2 




1 


l.OOOOe+00 


l.OOOOe+00 


l.OOOOe+00 


l.OOOOe+00 


2 


9.8341e-01 


1.0141e+00 


1.0188e+00 


1.0209e+00 


3 


1.9445e+00 


4.8168e+00 


1.2366e+02 


1.5746e+03 


4 


9.9933e-01 


9.9477e-01 


1.0340e+00 


1.0198e+00 


5 


9.8311e-01 


1.0596e+00 


1.0502e+00 


1.0290e+00 


6 


1.9989e+00 


4.9811e+00 


1.2663e+02 


1.6053e+03 






a = 


.25 




1 


l.OOOOe+00 


l.OOOOe+00 


l.OOOOe+00 


l.OOOOe+00 


2 


9.7862e-01 


9.6075e-01 


9.9969e-01 


9.7553e-01 


3 


1.0331e+00 


1.0580e+00 


3.3177e+00 


1.5299e+01 


4 


9.9450e-01 


9.7818e-01 


9.7776e-01 


9.8693e-01 


5 


9.5446e-01 


9.7849e-01 


1.0289e+00 


1.0236e+00 


6 


1.1033e+00 


1.1611e+00 


3.3920e+00 


1.5382e+01 



Table 11.1: Loop normalization constants. 



11.2 J = and the effect of torelons 

How well does this technique work? With only 1000 measurements, we are dealing with 
really low statistics. That these measurements are separated by 40 thermal sweeps each 
rather than the conventional 10 should help reduce the autocorrelations a bit, but we 
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11.2 J = and the effect of torelons 



really should not expect too much an exact agreement with earlier lattice results [1]. At 
this /3 and lattice size, 6 and 16 respectively, the mass of the 0^+ has been determined to 
be 1.193(18) lattice units, using a standard notation for the error ~ this is the number we 
are shooting for. Let's look at our correlation functions and see how we do. 

In Table 111.21 we have correlation functions before torelon subtraction. For lower a, 
the correlation function is overcome by noise around t=3, at least for the larger operators. 
At higher a we have a good correlation function for the size 4 operator - this operator 
seems to be our best for extracting the 0++ mass. The corresponding masses, and their 
errors, and shown in Table [TTT51 

After the torelon cancellation has been performed, we get the correlation functions in 
Table 111.41 and the masses in Table 111.51 We see that the correlation functions are better 
behaved by one more lattice spacing, lending more credibility to these mass estimates. In 
general, there is more consistency when the torelons are canceled. Also, the error on the 
mass estimates has been reduced. Looking at the second column of effective masses, we 
see that our best estimate of the 0"'"+ goes from 1.09(10) down to 1.06(08). We might 
even be seeing a mass plateau, as the next effective mass is 1.04(24), but saying that 
could be stretching our interpretation of the data by more than a little bit. While slightly 
low numerically, these estimates do agree well, within errors, with the canonical value of 
1.193(18). 
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t = 





1 


2 


3 


size 




a = 


.1 




3 


1 


2.5339e-01 


6.7977e-02 


1.3958e-02 


4 


1 


2.0777e-01 


4.8134e-02 


-3.2623e-03 


5 


1 


2.2440e-01 


5.3068e-02 


-3.4754e-03 


6 


1 


2.5338e-01 


7.1044e-02 


9.0107e-03 






a = 


.15 




3 


1 


2.5854e-01 


7.0588e-02 


1.4937e-02 


4 


1 


2.2336e-01 


5.3347e-02 


-2.1334C-03 


5 


1 


2.3454O-01 


5.5890C-02 


-2.3498C-03 


6 


1 


2..')8<)2o-()l 


7.3()78o-()2 


!).72()r)o-03 






a = 


.2 




3 


1 


2.6267e-01 


7.3012e-02 


1.5629e-02 


4 


1 


2.4336e-01 


6.2958e-02 


8.0794e-04 


5 


1 


2.4717C-01 


5.9903C-02 


-1.8301C-03 


6 


1 


2.6206e-01 


7.3653e-02 


9.2552e-03 






a = 


.25 




3 


1 


2.2035C-01 


5.9404C-02 


8.9081C-03 


4 


1 


2.4352e-01 


8.2043e-02 


1.9869e-02 


5 


1 


2.3038e-01 


6.4103e-02 


4.1644e-03 


6 


1 


2.2014e-01 


5.1320e-02 


-6.0157e-03 



Table 11.2: Correlation functions for 0++ before torelon removal. 



t = 


1 


2 


3 


size 


a = 


.1 




3 


1.37(03) 


1.32(10) 


1.58(50) 


4 


1.57(03) 


1.46(12) 




5 


1.49(03) 


1.44(13) 




6 


1.37(03) 


1.27(11) 


2.06(88) 




a = 


.15 




3 


1.35(03) 


1.30(10) 


1.55(45) 


4 


1.50(03) 


1.43(12) 




5 


1.45(03) 


1.43(13) 




6 


1.35(03) 


1.27(10) 


2.02(79) 




(\ = 


.2 




3 


1.34(03) 


1.28(10) 


1.54(45) 


4 


1.41(03) 


1.35(11) 




5 


1.40(03) 


1.42(13) 




6 


1.34(03) 


1.27(10) 


2.07(80) 




a = 


.25 




3 


1.51(04) 


1.31(14) 


1.90(86) 


4 


1.41(04) 


1.09(10) 


1.42(38) 


5 


1.47(04) 


1.28(14) 




6 


1.51(04) 


1.46(17) 





Table 11.3: Masses for 0"*^+ before torelon removal. 
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11.2 J = and the effect of torelons 



t = 





1 


2 


3 


size 




a = 


.1 




3 


1 


2.5339e-01 


6.7980e-02 


1.3961e-02 


4 


1 


2.2010e-01 


4.9200e-02 


8.4372e-03 


5 


1 


2.2544e-01 


5.4961e-02 


-2.2496e-03 


6 


1 


2.5348e-01 


7.1141e-02 


9.0913e-03 






a = 


.15 




3 


1 


2.5858e-01 


7.0627e-02 


1.4970e-02 


4 


1 


2.3390C-01 


5.8306C-02 


1.1263C-02 


5 


1 


2.3729O-01 


6.0708(;-02 


9.6583C-04 


6 


1 


2..')9.').")o-()l 


7.:-!749o-()2 


l.()272r-02 






a = 


.2 




3 


1 


2.6312e-01 


7.3452e-02 


1.6011e-02 


4 


1 


2.5080e-01 


7.0270e-02 


1.4823e-02 


5 


1 


2.5306C-01 


7.0020C-02 


6.0936C-03 


6 


1 


2.6515e-01 


7.7361e-02 


1.2315e-02 






a = 


.25 




3 


1 


2.2624C-01 


6.8414C-02 


2.0320C-02 


4 


1 


2.4196e-01 


8.3649e-02 


2.9644e-02 


5 


1 


2.3288e-01 


7.2597e-02 


1.7940e-02 


6 


1 


2.2597e-01 


6.5776e-02 


1.1193e-02 



Table 11.4: Correlation functions for 0++ after torelon removal. 



t = 


1 


2 


3 


size 


a = 


.1 




3 


1.37(03) 


1.32(10) 


1.58(50) 


4 


1.51(04) 


1.50(18) 


1.76(84) 


5 


1.49(03) 


1.41(12) 




6 


1.37(03) 


1.27(10) 


2.06(87) 




a = 


.15 




3 


1.35(03) 


1.30(10) 


1.55(45) 


4 


1.45(03) 


1.39(14) 


1.64(66) 


5 


1.44(03) 


1.36(10) 




6 


1.35(03) 


1.26(10) 


1.97(75) 




a = 


.2 




3 


1.34(03) 


1.28(10) 


1.52(40) 


4 


1.38(03) 


1.27(10) 


1.57(52) 


5 


1.37(03) 


1.28(08) 




6 


1.33(03) 


1.23(09) 


1.84(59) 




a = 


.25 




3 


1.49(04) 


1.20(07) 


1.21(32) 


4 


1.42(04) 


1.06(08) 


1.04(24) 


5 


1.46(04) 


1.17(08) 


1.40(41) 


6 


1.49(04) 


1.23(08) 


1.77(65) 



Table 11.5: Masses for 0+"'" after torelon removal. 
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11.3 J = 2: Coupling to higher spin states 



11.3 J = 2: Coupling to higher spin states 

Now that we are confident in our techniques for producing spin correlation functions 
and for removing the torelon contributions, we turn our attention to states with higher 
spin. The next hghtest state has 3=2. These operators have the phase contributions given 
by Equation 110.161 This operator is our first real test of our rotational phase method. 
For brevity, we will only present the results after torelon cancellation has been applied. 
The conventional estimate for the mass of the 2++, at these parameter values, from [1], is 
1.80(8). 

Table [TTTn] displays our correlation functions, and Table [TT771 the corresponding masses. 
These correlation functions fall to zero much faster, and so we can only go out to t=2 before 
statistical noise dominates. This state should be heavy enough that we can extract our 
estimate from the first column of effective masses. Also shown in the last column of the 
mass table is the overlap with the corresponding operator. These overlaps should 
tell us how well our operator is coupling to the 3=2. The lowest estimate of the mass is 
1.63(06); the next lowest is 1.68(05). These states have an overlap of about 1.5%. The 
third lowest estimate is 1.79(07), again with overlap of 1.5%. This amount of overlap is 
not sufficient to lower a mass of 1.8 down to 1.6: 

(1.5 X 1.2 + 98.5 X 1.8)/100 = 1.79 . (11.3) 

We might be wary of accepting the estimate coming from a = .25, as the only correlation 
function not hitting the statistical noise at t=2 is the size 4 operator giving us our lowest 
mass. Nonetheless, these are all in the 2a ballpark of 1.80(8), even if a little low. Thus, 
we do conclude that this operator is coupling effectively to the 3=2. Perhaps future 
calculations with better statistics will overcome this minor discrepancy. 
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11.3 J = 2: Coupling to higher spin states 



t = 





1 


2 


size 




a = 


.1 


3 


1 


8.0769e-02 


1.1711e-02 


4 


1 


9.3026e-02 


1.8498e-02 


5 


1 


8.0629e-02 


7.3778C-03 


6 


1 


9.8041e-02 


3.5999e-03 






a = 


.15 


3 


1 


8.4990C-02 


1.0426O-02 


4 


1 


1.2780e-01 


3.0552e-02 


5 


1 


9.8406e-02 


1.3467e-02 


6 


1 


1.0109e-01 


1.6209e-03 






a = 


.2 


3 


1 


8.6963e-02 


5.0222e-03 


4 


1 


1.8721C-01 


4.8227C-02 


5 


1 


1.3880e-01 


2.6834e-02 


6 


1 


9.9971e-02 


-4.0202e-03 






a = 


.25 


3 


1 


1.2540C-01 


-9.1018C-03 


4 


1 


1.9512C-01 


8.5575C-03 


5 


1 


1.6615C-01 


-2.0638O-05 


6 


1 


1.2756e-01 


-6.6520e-04 



Table 11.6: Correlation functions for 2++ after torelon removal. 



t = 


1 


2 


1 < 02\0o > 1 


size 


a ~ 


.1 




3 


2.52(08) 


1.93(05) 


.005(4) 


4 


2.37(07) 


1.62(03) 


.007(3) 


5 


2.52(09) 


2.39(07) 


.008(6) 


6 


2.32(06) 




.006(7) 




a = 


.15 




3 


2.47(08) 


2.10(50) 


.005(4) 


4 


2.06(06) 


1.43(19) 


.011(3) 


5 


2.32(08) 


1.99(42) 


.011(5) 


6 


2.29(06) 




.007(7) 




a = 


.2 




3 


2.44(08) 


2.85(1.15) 


.004(4) 


4 


1.68(05) 


1.36(13) 


.016(4) 


5 


1.97(08) 


1.64(23) 


.014(5) 


6 


2.30(06) 




.008(6) 




a = 


.25 




3 


2.08(07) 




.009(6) 


4 


1.63(06) 


3.13(1.11) 


.015(7) 


5 


1.79(07) 




.015(6) 


6 


2.06(07) 




.010(4) 



Table 11.7: Masses for 2++ after torelon removal. The last column is the overlap with the 
0++. 
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11.4 J = 4: The acid test 



11 A J = 4: The acid test 

Finally, we come to the operator we have all been waiting for: the J=4. Our hypothesis is 
that earlier lattice calculations have niisidentified the pseudoscalar as being J=0 when it 
in fact has J=4 (and hence is not a pseudoscalar). The state in question is the ^ with a 
mass of 2.10(33), from [1]. Our technique as developed so far only produces PC quantum 
numbers of H — h. However, gauge theory in D=2+l dimensions posses the curious property 
that in general states with J ^ come in degenerate parity doublets . Thus, while we 
might measure the mass of a state 4++, we may infer that there exists a state 4 ^ with 
equivalent mass. 

The correlation functions for J=4 are presented in Table 111.81 and the corresponding 
masses in Table [TTT^ Again, the overlap with the 0++ is shown in the last column of the 
mass table. This column is the one we are most interested in. Whatever operator has the 
lowest number here is our best estimate for the 4++. Scanning down these numbers, one 
state jumps out, the size 6 at a = .2. This state has an overlap of less than 1%! The next 
best operators have overlaps in the 5% range, too much for us to consider their masses 
as accurate. The mass of this state is 2.13(04). Let us emphasize that this error is the 
statistical error resulting from the jackknife analysis, and does not include the unknown 
systematic effects, which may be far greater. Nonetheless, it is pleasing to have such a 
tight statistical error resulting from a new technique of calculation. While this operator 
couples to the 4++ state, we invoke parity doubling to claim that there exists a state 
with quantum numbers 4 ^ at the same mass. Thus, we claim that the conventional 
assignment of ^ to the state measured with a mass of 2.10(33) is mistaken, and that 
this state really posses spin 4. Conventional lattice gauge theory needs to reevaluate how 
spin is assigned within the various symmetry channels in light of this result. 

These final results are shown in Figure [Tl. 21 
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11.4 J = 4: The acid test 



i 



X from Teper [1 8] 
o from calculation 



2 4 

J 



Figure 11.2: The results of the lattice calculation compared to the masses from Teper [1]. 



t = 





1 


2 


size 




a = 


.1 


3 


1 


8.4252e-02 


5.2959e-03 


4 


1 


1.1608C-01 


1.7138C-02 


5 


1 


6.4503e-02 


8.6513e-03 


6 


1 


9.8620e-02 


1.7948e-02 






a — 


.15 


3 


1 


6.5045e-02 


-2.7807e-03 


4 


1 


1.1552e-01 


1.7367e-02 


5 


1 


8.1668C-02 


1.7576C-02 


6 


1 


1.0721e-01 


2.0704e-02 






a = 


.2 


3 


1 


5.0137C-02 


-1.0371C-02 


4 


1 


1.3613e-01 


2.6607e-02 


5 


1 


1.3583e-01 


3.8166e-02 


6 


1 


1.1937e-01 


2.4707e-02 






a = 


.25 


3 


1 


6.9668e-02 


-1.7408e-02 


4 


1 


1.2066C-01 


7.5657C-03 


5 


1 


1.2103e-01 


1.7527e-02 


6 


1 


1.2769e-01 


4.7764e-05 



Table 11.8: Correlation functions for 4++ after torelon removal. 



2 

i 1.5 
a 

1 
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11.4 J = 4: The acid test 



t = 


1 2 


1 < 04\0o > 1 


size 


a = .l 


3 


2.47(08) 


.39 


4 


2.15(05) 1.91(28) 


.35 


5 


2.74(09) 2.01(99) 


.27 


6 


2.32(04) 1.70(53) 


.13 




a= .15 


3 


2.73(10) 


.27 


4 


2.16(06) 1.89(27) 


.24 


5 


2.51(07) 1.54(46) 


.17 


6 


2.23(04) 1.64(47) 


.05 




a = .2 


3 


2.99(12) 


.11 


4 


1.99(05) 1.63(21) 


.11 


5 


2.00(05) 1.27(18) 


.14 


6 


2.13(04) 1.58(39) 


.007 




a = .25 


3 


2.66(08) 


.03 


4 


2.11(07) 2.77(93) 


.38 


5 


2.11(05) 1.93(28) 


.05 


6 


2.06(05) 


.49 



Table 11.9: Masses for 4++ 
0++. 



after torelon removal. The last column is 



the overlap with the 



Chapter 12 



Conclusions 



As wc have seen, pure gauge theory does possess a rich striicture, even without the com- 
pUcation of quarks. Understanding the properties of the pure gauge field is vital if we are 
going to reach a deeper understanding of our universe. These gauge fields can produce 
quite interesting bound states: glueballs, gluelumps, and even glueknots. Even when re- 
stricting our attention to D=2-|-l, these glueballs in Flatland can intrigue and surprise us. 
Reaching towards a better comprehension of glueballs gives us insight into the structure 
of both mesons and baryons. 

The original model of Isgur and Paton was quite successful in explaining meson struc- 
ture. Its application to glueballs predated any experimental evidence for such states. 
Nonetheless, its rough predictions agreed well with the lattice data available at the time. 
While the extent of lattice calculations has improved greatly, not much work had been 
done to refine the model. Indeed, the most interesting refinements, such as the inclusion of 
a dissipating string tension and the possibility of the k-string, have required the input of 
research conducted within the last two years. Now, with a wealth of lattice data spanning 
several values of N from 2 to oo, in both 2 and 3 dimensions, we can test the predictions 
of these extensions against the lattice data. 

Several refinements to the original model have been proposed. The inclusion of an 
effective elasticity accords with the evidence coming from other forms of flux tubes and 
promotes the Liischer correction term from a constant to a free parameter for closed flux 
loops. Various mechanisms to produce states with negative charge conjugate eigenvalue 
have been developed, from a simple direct mechanism of fiux orientation reversal to more 
complicated mechanisms requiring the existence of higher topological contributions or of 
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more exotic types of gauge field flux tube. The properties of this flux tube at very smafl 
distances, where the idea of an infinitesimaUy thin flux tube must break down, has been 
incorporated. While still non-relativistic in character, these models are quite successful at 
reproducing the structure of the lattice spectrum. The best of these models, developed in 
Chapter [21 incorporates all these features: the k-string model with dissipating flux tube 
and effective elasticity. It can match the lattice spectrum to well within ±5% across a 
range of spins and excitations. 

One feature of all these models is that they disagree spectacularly with one particular 
state, traditionally assigned quantum numbers ^. No amount of parameter manipula- 
tion can bring that state into agreement without ruining the agreement across the rest 
of the spectrum, simply because any flux tube model must have several costly, ie heavy, 
phonons to produce a spin cancellation while providing negative parity. However, they 
all make a prediction that there exists a state with quantum numbers 4 ^ at just around 
that state's mass. Hence we are led to the hypothesis: "The lattice state traditionally 
labelled by ^ really possesses a spin of 4." Such an hypothesis cannot be automatically 
discounted owing to the mod4 spin ambiguity of operators on a cubic lattice. 

To test this hypothesis, a new approach to constructing operators on a cubic lattice 
is developed, starting in Chapter (5] This approach borrows from techniques of Greens 
function inversion already used in fermionic calculations, but applies the inversion to the 
full lattice, not just to all paths emanating from one particular site. The munerical cost of 
this technique is more expensive than traditional methods, but the potential to overcome 
cubic group spin ambiguities makes its pursuit worthwhile. As this technique receives 
further exploration, methods to tame the numerics will be developed, as in Section [10.21 

The result: there does exist a state with quantum numbers 4 ^ precisely at the mass 
which the lattice community calls the ^. This is a beautiful result, conflrming by a 
lattice calculation an unconventional hypothesis coming from a dynamical model of flux 
tube glueballs. The agreement is remarkable, comparing a mass of 2.13(04) with a mass 
of 2.10(33) respectively. At the end of the day, such close agreement strengthens both 
aspects of the hypothesis, giving much credibility to the flux tube model of glueballs as 
well as inspiring further research into the method of Greens function inversion. 



Appendix A 

Bayesian Methods of 
Inferential Statistics 



When comparing a model to a set of data points, the most appropriate method to use is 
that of inferential statistics [38,39]. This method provides means to estimate the values 
of the parameters that are most likely, given data points to match, as well as to compare 
the relative likelihood of several competing models. What we are given is a set of data 
{D} of lattice measurements from which we are to infer the most probable values {a^} 
of the parameters of the underlying theory. What we actually compute is the probabil- 
ity distribution function as a function in multi-dimensional parameter space. The joint 
distribution of X and Y is given by Bayes Theorem as 

Mx|r,/) = 5ffli^, ,A.i, 

where X and Y are propositions conditioned on the background knowledge /. When these 
variables represent our hypothesis and data, we have 

pr {hypothesis] data, I) oc pr{data\ hypothesis, I) x pr {hypothesis] I). (A. 2) 

The Ihs is called the "posterior" and is the function we wish to maximize. The first factor 
on the rhs is called the "likelihood" of the data given the hypothesis (and /). The second 
factor is the "prior" , and is often a source of confusion for beginning Bayesians. Physicists 
often assign a uniform prior (sometimes when not appropriate!), reducing the problem 
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to the familiar one of maximum likelihood parameter estimation. However, the Principle 
of Insufficient Reason, as stated by Bernoulli, can prescribe more useful priors in most 
circumstances. The omitted factor, the "evidence", is only a normalizing constant when 
pcforming parameter estimation, evaluated by normalizing the posterior, but reappears 
when we become interested in model selection. 

For uniform priors, the method reduces to the familiar one of maximum likelihood 
estimation [34]. The likelihood, under certain assumptions, can be taken to be x^, a 
nonnegative quantity. As the probability is defined as the exponential of x^j we often 
work with the log of the posterior, which would simply be our likelihood function. 

For future reference, we state the means by which we remove unwanted propositions 
from a joint distribution. Continuous marginalization is given by 

/ + 00 
pr{X,Y\I)dY (A.3) 
-oo 

which goes over in the trivially ultimate discrete case to 

pr{X\T) = pr{X,Y\T) + pr{X,Y\I). (A.4) 

When the posterior is sharply peaked around a unique maximum, we may make a local 
quadratic approximation to determine the uncertainty in a fitted parameter by fitting 
a local Gaussian to the peak. This case corresponds to the usual notion of confidence 
intervals. When the posterior is not sharply peaked, or even multinodal, such a gross 
approximation is not adequate. Even so, in this case we may still define Bayesian intervals 
to represent our state of knowledge of the parameters given the data at hand. 



Appendix B 



Review of the 



Isgur-Paton 



Model in Two 



Spatial 



Dimensions 



While perturbative QCD works demonstrably well for high-momentum scattering experi- 
ments, the phenomenology of hadrons lies beyond perturbative methods. Quark models, 
on the other hand, give a reasonably good accounting of the hadronic spectrum. Pull 
QCD, though, must have a spectrum more complicated than the quark model, as the self- 
coupling of the gluonic sector predicts bound states devoid of any quarks. These states, 
called "glueballs" or "gluelumps" , can be inferred from the spectrum of the pure Yang-Mills 
theory corresponding to quarkless ("quenched") QCD. Two avenues for approach present 
themselves: weak-coupling models or strong-coupling models. Weak-coupling models are 
epitomized by the MIT-bag model [16] , where "constituent gluons" are postulated, in anal- 
ogy with the familiar "constituent quarks" , to live within a confining potential bag. While 
these models can reproduce some features of the Yang-Mills spectrum, many issues are 
unresolved. Additionally, these models actually violate the "no-hair" theorem [17], at least 
at the semi-classical level. Turning to strong-coupling models, Isgur and Paton developed 
a model [6] where the gluonic degrees of freedom condense into tubes of chromoelectric 
flux between static quarks. The quarkless extension of this point of view envisages pure 
glue states as a closed loop of flux-imagine annihilating the q and g of a meson leaving 
the flux tube behind as ring of flux. Extensions to this naive flux tube model are possible. 
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with higher representation flux and more compUcated topologies allowed, as long as fliix 
junctions are properly treated and final states remain gauge invariant. 

Without the presence of the static quarks to provide a natural adiabatic limit, Isgur 
and Paton supposed that the low-lying spectrum will be described by the simplest state: 
a closed loop of fundamental flux. A detailed application of the Isgur-Paton model to the 
case of D=2-|-l and SU(2) was performed in [28]. The flux tube is described by polar 
coordinates (r, 6). For small fluctuations about a circle of radius tq, we expand the radial 
variable in Fourier modes which we identify as phonons. These phonons contribute to 
the potential energy of the Schroedinger equation for the radial flux tube when we apply 
an ad hoc Born-Oppenheimer approximation to absorb the "fast" phonon modes. The 
resulting Hamiltonian is diagonalized to give the spectrum. The traditional Isgur-Paton 
model contains a cutoff parameter /, as well as the string tension a which is set equal to 
unity. 

For SU(N>3) the flux tube has a direction, so a circular string will have one of two 
orientations. These we shall label by L and R. In 2 space dimensions one cannot change 
the orientation of a loop by a rotation so there is no link between the rotational and 
charge conjugation properties of a state (in contrast to the case in 3 space dimensions). 
Moreover the model provides no means by which an L loop will transform into an R loop. 
Thus the states which are labelled by the orientation, O = L ov R, of the loop, by the 
phonon occupation numbers, {n^}, and by the radial quantum number, n^, are energy 
eigenstates in the model and form a complete set. We may write them as \0; {n^}; n^). 
Since the theory is invariant under charge conjugation and parity we can also label the 
energy eigenstates by these more physical quantum numbers. What are these states? 

The action of charge conjugation, C, is to reverse the orientation of the glue loop 
without changing its modes in any way: 

l^^inrnhnR) \R;{n^};nR). (B.l) 

Moreover for all quantum numbers a state and its charge conjugate are different. Thus 
when we apply 1 ± C to one of our energy eigenstates we get a pair of non-null degenerate 
states of opposite quantum numbers. Hence the masses in the C = + and C = — sectors 
are identical. It is clear that this degeneracy simply follows from the fact that there is no 
mixing in the model between loops of opposite orientation. 
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Figure B.l: Transformation mider parity and charge conjugaion. 

Turning now to parity, P, we recall that in 2 space dimensions it transforms (x, y) 
(x, —y) (up to a rotation) and so changes the sign of angular momentum. Thus it will 
reverse the helicity of each phonon. It is also clear that it reverses the orientation of a 
closed flux loop. Thus 

|L;{n±};n;j) ^ (B.2) 

These two states are different but degenerate. Thus when we apply 1 ± to form the two 
parity eigenstates we get a pair of non-null degenerate states. That is to say, the mass 
spectra in the P = + and P = — sectors are identical. For J ^ this degeneracy is in 
fact quite general, in 2 space dimensions. However the general argument breaks down for 
J = and here it is the special features of the model that lead to such a degeneracy. 
A pictorial representation of the action of these transformations is given in Figure lB.ll 
From the above we see that an energy eigenstate definite P and C can be written (up 
to a normalisation) as 

^P^c - i\L; ni) + rjp\R; n^)) + j^^m ^t) + Vp\L; n^)) (B.3) 

where r/p — ±1 is the eigenvalue for parity and 7]c — ±1 is the eigenvalue for charge 
conjugation. (We have abbreviated our notation for the states in an obvious way - in 
particular by suppressing the radial quantum number, np, which is unchanged by both P 
and C.) For example, the lowest lying J — 2 state has a single m — 2 phonon, and so the 
P = +,C = —,J = 2 ground state is 

i;j=2 = {\L; n+ = l) + \R; = 1)) - {\R; n+ = 1) + |L; = 1)) (B.4) 

From Equation IB. 31 we can readily infer some qualitative features of the SU(A^ >3) 
mass spectra: 
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1) The very lightest states are those with no phonons at all (and with riji = 0). There are 
two such states, the 0++ and the , and they will be degenerate. 

2) Since each of the component states in Equation IB. 31 is itself an energy eigenstate and 
they are all degenerate, this implies that for a given J the four sectors labelled by the 
four combinations of P, C are degenerate except where one of the linear combinations 
is null. The latter can only occur if the phonon content is identical for both helicities, 
{"-ml — {"-ml: ^his Can only occur for J = 0. So, for example, if we take the state 
with no phonons at all and try to form a ^ state degenerate with the the 0++ ground 
state we find, from Equation IB. 3[ 

V^o-^ = i\L; 0) - \R; 0)) + (|i?; 0) - \L; 0» = (B.5) 

It is easy to see that the lightest ^ arises for — l,n2 — 2 and so is a (very) heavy 
state. 

All the above is for SU(A^ >3). The SU(2) spectrum has no C — sector which, in the 
string model, is embodied in the fact that there is no orientation on the flux loop. So if we 
identify L and R in the above equations we obtain expressions for the energy eigenstates. 
Once again we have parity doubling for J ^ and for that part of the J = spectrum that 
contains differing positive and negative phonon contents. The following table summarizes 
the assignment of spin and parity quantum number available in D=2+l for the Isgure- 
Paton model. 



M 







0+ 


2 


2± 


3 


3± 


4 


0+,4± 


5 


1±,5± 


6 


0+,2±,6± 


7 


l+,3±,7± 


8 


0±,2±,4±,8± 


9 


1±,3±,5±,9± 




etc. 



Table B.l: Allowable quantum numbers for each value of M. 



Appendix C 

Programs Used for the 
Calculations 

All calculations were performed using the Matlab programming environment. Some 
advantages of Matlab include the ease with which complex algorithms may be altered 
and its more mathematical programming environment. The inherent matrix data type 
lends itself well to prototyping logical algorithms without having to change large swaths 
of dedicated C or Fortran. Also, the mathematics behind the calculation are more easily 
seen using Matlab programs. The crucial disadvantage is speed. While Matlab does 
run remarkably fast, its inherent overhead slows down production calculation. It simply 
cannot compete with highly optimized and parallelized code. For the purposes of these 
calculations, it was the appropriate tool to use, but for any future work (on the lattice 
calculations), optimized C is definitely the way forward. Note, the SuperLU code [54] was 
modified slightly to create a .mex file which handles complex values. 
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% drives minimization of doglueStoo et al 

% here IP model with f , g compared to C=+ sector for all N 



% modified Sept 3, 1998 to reflect 

global data dataerr 



■ SU(infinitY) data 



data ( 
dataerr ( 
data ( 
dataerr ( 
data ( : 
dataerr ( 
data ( : 

data ( : 

data (1 
data (1 

%datae 



2, 1) = [4 . 718 7 . 
1 : 2, 1 j = [ . 043 . 1' 
,2)-[4.329 7. 
1:2,2)-[.041 .] 
,3t-[4,236 7, 
l:2,3t-[.050 .] 
,4)-[4.184 7, 
l:2,4)-[.055 
,5)-[4.065 6. 
l:2,51-[.055 
1) -8 .15; data (1, : 
4) -7 .85; data (1, ; 



2 10.42;6.83 

. 34; . 10 inf 

3 10.22;6.52 
.24; ,09 

5 9.91;6.38 
,36;. 13 .20 

9 10.26;6,20 
.50;. 13 -18 

8 9.98;6.18 
.25; . 13 

2) =8 .23; data (1, 3, 

5) =7 . 99; 



inf]', 

8,59 

inf]', 



1) 15; dataerr (1,3, 2)=. 15; dataerr (1, 3, 3) - 
4)-. 22; dataerr[l,3,5j=.22; 
■r (1, 1, :) -inf .'ones (1,1,5);; 
■=sqrt {dataerr. ^2 + ( . 05 . -data) . ~2) ; 



dataerr ( 
data ( : 
dataerr ( 
data ( : 
dataerr ( 
idataerr 



, 1) = [ 9 . 95 
1:2, 1) =[ . 32 
, 2 j = [ 9 . 30 
1 : 2, 2 j = [ . 25 
,3j-[9.31 
1:2, 3)-[.28 
, 4t-[9,19 
1:2, 4t-[.29 
,5)-[9,02 
l:2,5t-E.30 
, 1, :)-inf .*oi 



dataerr^sqrt (datat 



.86 11.13;0 
14 .42; inf inf 
.36 10.19;G 
11 .27; inf 

.86 10.85;0 
18 .55;inf 
,18 10,28;0 
16 .34;inf 
,89 10,06;0 
21 .40;inf 
(1,1,5);; 
(,05, -data) ."2); 



data ( : 

data ( : , 
dataerr ( 
data( 



2. 1) -[0 
,1:2, 1) = [inf 

2.2) -[6,48 
,l:2,2)-[.09 

2, 31=[6.271 
, 1 : 2, 3) = [ . 095 
2, 4) = [ 6 . 03 
, 1 ; 2, 4) = [ . 13 



chisql=doglueIPnewk2 {f gaout , Ms, N {dk) 
M - [0 2 3 5 4]; 

[chisq2, massout , mc, oP, oM] =doglueIPnewk2 (f gaout 
(dk) ] ,1, 1) 

fout- [' b5newk2N' nura2str (dk) ' .mat' ] ; 
eval ( [' save ' foot ] ) 
end % for dk 

lerr^lasterr 

fout=['b5newk2N' nijm2str (dk) ' err .mat' ] ; 
eval ( [ ' save ' f out ] } 



function [chii 



% calculates masses for glueballs 
% with C-splitting 

by Robert W. Johnson 



Lrgin<3, sig-2; end 



3-ga (1) ; 

^ga(2) ; 
b-ga(3) ; 



,Ms, ",datal,dataerrl) 
% chisql^dogluelPadj (f gaout, Ms, ' ' , datal, dataerr 1 ) 

[ f gaout, opt] =fmins (' doglueIPnewk2' , f gain, opt, [ ] , Ms, N (dk) , 



datal, dataerrl) 



l,N{dkJ, ['b5newk2N' num2.=;rr 



lap, olapm] ^dogluelPnewk (ga, Ms, N, f name, data, dat ai 
Tsgur-Paton model 



Lrgin<4, fname - 



argin<5, data = 
argin<6, datae: 

argin<8, neig 



(12, 3) ; end 
nes (12, 3) ;end 
1000; end 



ros(3*n,neig) 
,neig) ; V=U; 

ig) ; 



%iters - 
P - zero; 
D = zero; 
HI = zeri 
H2 - zeri 



lloc(3*n, 3*n,21'nj ; 
eros (neig, 1) ; 
ros (neig, 1) ; 
OS (neig, 1) ; 
zeros (neig, 1) ; 
OS (neig, neig) ; 
OS (neig, neig) ; 
ros (neig, neig) ; 
ros (neig, neig) ; 
ros (neig, 1) ; 
(k-l)*2+l:k*2; col ■ 



dataerr (1, j, 2) -.2 6; aatacrr ( 1 , 3, 3 ) = . 30 ; 
dataerr (1,3, 4)=. 41; dataerr (1,3,5)=. 65; 
%dataerr (1, 1, : ) =inf . *ones ( 1 , 1 , 5) ; ; 
dataerr^sqrt (dataerr . '2 + ( . 05. *data) . "2 



boff = b'( 



, 1) ; 



data ( : 
data ( : 



dataerr ( 
data( 
dataerr ( 
%datae 



2, 1)=[0 
1:2, 1) = [inf 
, 2) = [10 . 52 
1:2,2)=[ .28 
, 3)-[10.35 
1:2, 3)-[.S0 
, 4)-[9.43 
1:2, 4)-[.75 
,5)-[9,47 
1:2,5)-[1.16 
: ) -inf . *ont 



0;0 
inf 
I 10.51 



.45; inf 
9. 93;0 
. 44; inf 
9.43;0 ' 



. 05 . *data) . "2) ; 



dataerr-sqrt (datai 

data- [dp; dl ; data] ; 

dataerr- [dperr; die; dataerr ] ; 

%data-data( [1 72 11 583 12 594 10],:,;); 
%dataerr-dataerr[ [1 7 2 11 5 8 3 12 6 9 4 10 ] , 
clear dl die dp dperr 

>I-[2 3 4 5 1/eps]; 

% begin program 

3-[.5:.5:l,5]; 
3-[-fliplr(g) g ]; 
a-[0: .25:2] ; 
4b-[0: .5:5] ; 
[G,A] -ndgrid (g, a) ; 



[Sip, Sk, rh, T] - maketrixIPk (n, M, g, W) ; 
Sm — [Sip 4*pi*c"2*spdiags (aof f , 0, n, n) 
spalloc (n, n, n) ; ... 

4*pi*c"2*spdiags (aof f , 0, n, n) Sip 
(boff,0,n,n) ; ... 

4*pi*c^2*spdiags (bof f , 0, n, n) spalloc 
[aoff,0,n,n) ; . , , 

spalloc (n, n, n) 4*pi*c'^2*spdiags (boff, 
n,n) Sk]; 

opt .tol=le-10; opt .di3p=0; 
[Dip, Dip] =eigs (Sip, neig, 0, opt ) ; 
[Uk, Dk] =eigs (Sk, neig, 0, opt ) ; 
[C-, D-] =eiqs { Sm, 2*neig, 0, opt ) ; 



Rip = Sip^Uip -Uip'Dip; 
Rk- Sk*Uk -Uk*Dk; 
Rm— Sm*Um -Um*Dm; 

U - [Uip Uk] ; 

V-U,*repmat{T, [1 2*neig] ) ; 

ma - l./sqrt(-trapz(rh, V.''2))' 

V - V*diag(ma) ; 

Vm=Um. *repmat (T, [4 2*neig] ) ; 
ma = 1 . /sqrt (-trapz (rh, Vm[l:n, 
'm (2 'n+1 : 3*n, : ) . '2) ... 

- trapz (rh, Vm(3*n+1 : 4*n, :). " 
Vm = Vm^diag (ma) ; 



4 *pi*c"2*spdiags 
spalloc (n,n,n) 



4*pi*c''2- 
4*pi*c''2 
2*spdiags 



.,n) 

.*spdiags 

ipdiags 
(aof f , 0, 



: ) . "2) - trapz (rh, Vm (n+1 : 2*n, : ) , -2) 



olap = U' *U; 
olapm — Urn' 'Um; 

% need to sort Um, Dm 
for kk— l:2*neig 
p (kk) — prod (ni 



(bandmax (Um (1 : 2*n, kk) ' , 1, 1) } ) 



s ( [length (g) , length (a) ] ) ; 



for dk-2:5 % N-2, 3, 4, 5, inf 
Ms-[0 2]; 



datal-squeeze (data { : , ; , dk) ) ; 
dataerrl-squeeze (dataerr { : , ; , dk) ) ; 



opt{l)-0; opt(3)-le-4; opt(2)-le-4; opt{14)-300; 

for ak— 1 : length (a) 

for gk=l ; length (g) 

for bk=l : length (b) 

chisq (gk, ak) — doglueIPnewk2 ( [g (gk) 



(ak) ] ,Ms,N (dk) , 



, datal, dataer 



zin-f ind (chisq ( : ) ~min (chisq (:))); 
[zg, za] -ind2 sub (size (chisq) , zin) ; 
fgain-[g(zg) a(za)]; 

[ f gannt, opt , 1 am, he.q s ] -con str (' dngi ueTP^d j' , f grii n, opt, [-1.5 0] , [ 1 . 5 



mm-fi-d (-P) ; 
V=Vm(:,m:^) ; 

massU=diag (Dm (mm, mm) ) ' ; 

Um{;,mm)-[]; 

Vm(;,mm)-[); 

Dm(:,mm)-[]; 

Dm(mm,:)-[]; 

massU - massU/(4*pi*c'-2) , 
massM — diag (Dm) ' / ( 4*pi*i 



massout ( ; , : , k) - [massU;massP;raassM] ; 

if -isemptyffname) , 

f nam— [ ' tmpM' int2str(M) char(fname) '.mat']; 

eval ( [ ' save ' f nam] ) ; 
end % if 

id ^ for k 



Figure C.l: Programs used for the extended Isgur-Paton model. 
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% need to sort masses to make chisq 
k currently fitting M=0,2 

mc=zero3(12,3); 

mc(l, l;3)^a330Ut (3,1:3,1); 

mc(2,l:2)-massout(3,l:2,2); 
% mc(3, :)-inassout (3,1:2,3) ; 
% mc (4, : )=massout (3, 1:2, 4) ; 

inc(5,l:2)=massout(3,l:2,2); 
% mc (6, : ]=massout (3, 1:2, 3) ,■ 
mc (7, 1 : 3 ) =massout (1 , 1 : 3, 1) ; 

mc (8, 1 : 2) =massoiit (1,1:2,2) ; 
% mc [9, : )-massout (1,1:2,3); 
% mcllO, :)-massout(l,l:2,4); 

mc(ll,l:2)-massout{l,l:2,2) ; 
% mc{12, :)-maS30ut (1,1:2,3); 

%remQves J=l from fit 
mc( [3 4 6 9 10 12] ,:) = [] ; 
data( [3 4 6 9 10 12] ,:)-[] ; 
dataerr([3 4 6 9 10 12], :)-[]; 

chisq - sum( { (data{:)-mc{ :) ) ./dataerr (: ) ) .^2) ; 

disp (lasterrj 
save chierr .mat 




ha- z - aa; 

% need values of spectrum at points around 
% mgl and mah lower and higher for g and a 
eval ( [ ' [z, mass] =doglueTPnewk2' . , . 

' ([gg-hg aa] , [0 2],n);']) 
mgl-zeros (size (dataerrl) ) ; 
mgl (1, :)-mass (3, 1:3,1) ; 
mgl (2, :)-mass (3, 1:3,2) ; 
mgl (5, :)-mass (3, 1:3,2) ; 
mgl (7, :)-mass (1,1:3,1); 
mgl (8, :)-mass (1,1:3, 2) ; 
mgl (11, : )=mass(l, 1:3,2) ; 
% removes J=l from fit 
mgl ( [3 4 6 9 10 12] ,:)-[] ; 

eval ( [ ' [z, mass] =doglueIPnewk2' , . . 

' ([gg+hg aa] , [0 2],n);']) 
mgh-zeros (size (dataerrl) ) ; 
mgh(l, : )-mass (3,1:3,1); 
mgh(2, : ) -mass (3,1:3,2); 
mgh(5, : ) =mass (3,1:3, 2) ; 
mgh(7, : ) =mass (1, 1 : 3, 1) ; 
mgh(8, : ) -mass (1, 1 : 3, 2) ; 
mghdl, :) -mass (1, 1:3,2) ; 
%removes J-1 from fit 
mgh( [3 4 6 9 10 12] ,:)=[] ; 

eval ( [ ' [z, mass] — doglueIPnewk2 ' ... 

' ([gg aa-ha], [0 2],n);']) 
mal-zeros (size (dataerrl) ) ; 
mal (1, :)-mass (3, 1:3,1) ; 
mal (2, : )-mass (3,1:3,2); 
mal (5, : )-mass (3,1:3, 2) ; 
mal (7, : ) -mass (1,1:3,1); 
mal (8, : ) -mass (1,1:3,2); 
mal (11, : ) -mass (1,1:3,2); 
"ijremoves J-1 frcm fit 
mal { [3 4 6 <• IC 12] , ; ) - " ] ; 

eval ( [ ' [ z, mass] — doglueIPnewk2 ' ... 

' ([gg aa+ha], [0 2],n);']) 
mah-zeros (size (dataerrl) ) ; 
mah(l, :)-mass (3,1:3, 1) ; 
mah (2, : ) -mass (3,1:3,2); 
mah (5, : ) -mass (3,1:3,2); 
mah (7, : ) -mass (1,1:3,1); 
mah ( 8 , : ) -mass (1,1:3,2); 
mah (11, : )-mass (1,1:3,2); 
% removes J— 1 from fit 
mah( [3 4 6 9 10 12] ,:)-[]; 

dg - (mgh - mgl) . / (2*hg) ; 
da = (mah - mal) . / (2*ha) ; 

derr-derr ( : ) ; 
dg-dg { : ) ; 
da-da ( : ) ; 

% elements of curvature matrix 

cur- [sum (dq . ^ dq . /derr . "2 ) sum (dg. 'da . /derr . '2) ; ... 




function [Sip, Sk, rh, T] =maketrixIPk (n,M, g,N) ; 

% n is discretization scale, M is phonon parameter 

% written by Robert W. Johnson 

persistent r c f 

Sb = spalloc [n,n, 5-n) ; 
k = spalloc (n, n, n) ; 
Fl - spalloc {n,n,nj; 
Fy - spalloc [n,n, n) ; 
kspeye = speye (size (kj ) ; 
3d = zeros (1, n) ; 
3dd= zeros {1, n) ; 

sigK=3in(2^pi/H) /sin(pi/N) ; 

li=l/(n + l); % interval length 

% create matrix 

rh - ( [n+1) ./j - 1) ; 

4 load r and c and f 
load rvhro 

%t=j.*h; 






j - ones(n,l); 

Sb = spdiagsdj -16'j SCj -16-*] j], -2;2, n, n) ; 
Sb(l,l)=31; Sb(n,n) - 29; % boundary condition 

Sd - spdiags[2*pi.*(l-exp(-f.*r)) .-2.*r,0,n,n) + (M + g) * ... 

spdiags (l./r,0,n,n) ; 
Sd=4*pi*c"2.-'Sd; 

Sk= spdiags (2*pi'-sigK. * (1-exp (-f. *r) ). ~2 . 'r, 0,n,n) + (M + sigK'g) * ... 

spdiags (1 . /r, 0, n,n) ; 
Sk-4'-pi'-c'2.*Sk; 

Sip = h"2/12-*k*Sb*k+Sd ; 
Sk - h^2/12/sigK*k*Sb*k+Sk; 




"J; clear unnecessary variables 

clear mah mal mgh mgl h z gg aa mass 

eval ( [ ' save ' fin ' err .mat' ] ) 

%end 




n-1000; 
li-l/ (n+1) ; 
j-l;n; 

t-j.*h; 

rho=I./t -I; 

r=zeros (size (rho) ) ; 

r(l)=fzero('rhof',rho[l)/2, [ ] , 0, rho (1) , f2, c) ; 
for z-2:n 

r (z)-fzero['rhof' ,r (z-1) , [] ,0,rho{z) ,f2,c) ; 


% computes errors 

% modified by hand to deal with TPver 

%for t=l;2 

for b=[5 0] 
% if t==3 
% ns=2:4; 
% else 
% ns-1; 

for n-3 

fin=['b' num23tr [bt ' newk2N' ... 

num2str (nj ] ; 
eval ( [ ' load ' fin] ) 

gg=f gaout ( 1 ) ; 
aa=f gaout ( 2 ) ; 
derr=dataerrl ; 
derr ( [3 4 6 9 10 12] ,:)=[] ; 

% select steps ize for differences 

h - eps-^ (1/3) /2; 

h - .01; 

z - gg + h; 

hg- z - gg; 

h - .05; 

7 - sfl + h: 











Figure C.2: Programs used for the extended Isgur-Paton model. 
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% Script dolatJ.m 

laster lattice simulation script 
k Robert W, Johnson 

% Global Variables 

ibal Psi Parameters Xc Xs Xp Xm 

'arameters 

:ould be read in from a file 



'D', 3, 
' Beta' , 
' LatSize 
'Start' , 
' StartSweep' , 
' Stagger' , ; 
' NumSweeps' , 
' Thermalize' 
' Sweeps' , 11 
'Bins', 20, 
'NumCM', 1, 
'NumOR', 5, 

'GlobalGT' , 
'CalcPlaq' , 
'ValPlaq' , ( 
'Accept' , 0, 
'LatRng' , 0, 
'SngRng', 0, 
'PlaqReal', 0, ,. 
' Plaqimg' ,0, 
' InFile' , ' therm' , 



% SU(Nj 

% D = Euclidean dimension 
% beta = 2N /g0^2 
% [X y t] 

% (0,l,2}-{load,cold,hot} 
% or load with configuration 
% sweep staggered update2s 
... % compound f 



.atii 



'OutFile', 
) ; 



'/v, 



. ■ . % thermal; 

% sweeps between measurements 

% bins for error anal; 
% Cabibbo-Marinari swt 
% cvcr-rclax sweeps p. 
% reunitarizc c-vcry nth compoi 
% global gauge transformation 
% calculates plaquettes 
% validate plaquettes 
% l=always accept 
% lattice random generator 
% scalar random generator 
% check when loading 
% check when loading 
... % filename to load 

r/freedisk/r johnson/b6124therm' 



to throw away 



ery nth compoi 



% filename to : 



sp (Parameters) 

,gger=Para:neters . Stagger; 
L=Parameters . LatSize; 
3=Parameters .D; 
Sl^Parameters .K; 
i^Parameters .Bins; 
PN^Parameters .NumSweeps; 
iump^Parameters , Sweeps; 
|outfile=char {Parameters .OutFile) ; 

-zeros ( [D*N. ^2 L"D] ) ; 



Kp-repmat{Xc, [1 D] J ; 

Kp(:,l)-sub2ind([I. !■ !■] , mod (Xs ( : , 1) , L) +1, Xs ( : , 2J , Xs ( : , 3) f ; 

!<p ( : , 2 j =sub2ind [ [L L L] , Xs ( : , 1) , mod (Xs [ : , 2) , L) +1, Xs ( : , 3 ) 1 ; 

: , 3) =sub2ind ( [L L L] , Xs ( : , 1 ) , Xs ( : , 2 j , mod (Xs ( : , 3) , L) +1 ) ; 

■leprat (Xc, [1 Dl ) ; 



performs calculations on saved configs 

global Psi D L N Xc Xs Xp Xm Xops 

nfile^' /home/wytham/r johnson/myfree/therms/b6116therm' ; 
outf ile=' /home/wytham/ r johnson/myf ree/ops7/ops7b6116' ; 

5val(['load ' infile 'l.mat']) 
je calcl6err.dat h -ascii 



Hsweeps=4000; 

kip=4; 



R-stapleJ(sites (k, :) , mu) ; 

Psi ( (mu-1) *N'^2+[1:N~2] , sites (k, :)) -update J (sites (k, :) , mu, R) ; 



m (hmc. Parameters , GlobalGT) 



% global GT 



% V now SU{N) 



V=rand ( [N"2 L''D1 ) ; 
nv=sqrt (sum(V.'-2) ) ; 
V=V./repmat(nv, [N^2 1]); 
for mu=l:D 

Psi ( (mu-1) *N-^2+[l :N~2] , : ) -ml (V,m2 (Psi ( (mu-l) *N~2+ [1 :N^2] , : ) , 
V(:,Xp(:,mu)))); 



for n=l : Parameters .NumOR 



% heatbath + OR 



ela 



for mu-l:D 

for k=l: Stagger 

R-stapleJ(sites (k, : ) , mu) ; 

Psi ( (mu-1) *N"2+[l;N"2j , sites (k, ;)) -ovrlaxj (sites (k, : ) , 

end 

end 

end % OR sweeps 

if -rem(hmc. Parameters .Reunit) % reun 

for k-l:D 

nPsi-sqrt [sum(P3i( (k-1) *N~2+ [1 :N''2 ] , :) ."2J ) ; 
Psi ( (k-lj *N''2+[1:N'2] , : )-Psi ( (k-l) *N"2+ [ 1 ;N''2] , ; ) . / ... 
repmat [nPsi, [N''2 1]); 

end 



% di3p{hmc) 
eval(['save ■ 
end % if 

tnd % for hmc 



tfile num2st r [hmc/dump j ' 



k-0:Stagger-l 

lites (k+1, :) -find (mod (sum (Xs' ) , Stagger) -=k) ; 



% setup gauge fields Psi 
litch Parameters . Start 
% loads configuration 

eval(['load ' char (Parameters . OutFile) ] ) 



% cold start 

Psi=zeros ( [D-*N. '2 L"D] ) ; 
Psi (N-'2 . * [1 :D] , : ) -1; 



Psi=rand(sPsi) ; 
for k-l:D 

nPsi-sqrt {sum(Psi ( (k-l) *N-2+ [ 1 :N''2] , :) .^2) ) ; 
Psi { (k-l) *K~2+ [1:N''2] , : ) -Psi ( (k-l) *N^2+ [1 :N'2] , : ) . / ... 
repmat (nPsi, [N''2 1]); 

=nd % switch 

for hmc—1 rParameters . Thermalize % thermalization sweeps 
for mu-l:D 

for k=l: Stagger 

R-stapleJ{sites{k, :) , mu) ; 

Psi{ (mu-1) *N^2+[1:N'^2] , sites (k, : ) ) -updateJ (sites (k, :) , 

end 

end 

for n=l : Parameters . NumOR % over-r 

for mu=l:D 

for k=l : Stagger 

R=stapleJ(sites[k, :), mu) ; 
Psi ( (mu-1) *N"2+ [1 :N"2] , 

end 



la=length(alph) ; 

set up operators 
.-20; 
L-16; 
D-3; 
H-2; 

ops- [2 3 4 5 6]; 
lop3-length(ops) ; 

^c-[l:L~3]'; 

Xs-repmat (Xc, [1 3] ) ; 
[Xs ( : ,1) ,Xs(:,2) ,Xs (: 

;<p=repmat (Xc, [1 D] ) ; 

:<P ( 
:<P ( 

!<p [ 



3) ]-ind2sub( [L L L],Xc); 



2ind{[L L L] , mod (Xs ( : , 1 ) , L) +1 , Xs ( : , 2 ) , Xs ( : , 3) ) ; 
, 2)=sub2ind( [L L L] , Xs ( : , 1 ) , mod (Xs ( : , 2 ) , L) +1, Xs ( : , 3) ) ; 
, 3)-sub2ind( [L L L], Xs ( : , 1 ) , Xs [ : , 2) ,mod [Xs [ : , 3) , L) +1) ; 



Km-repmat (Xc, [1 D] ) ; 

>tm [ : , 1) -sub2ind ( [L L L] , 

>;m ( : , 2] =sub2ind { [L L L] , 



( [- 



mod(Xs (:, l)-2,L)+l,Xs (:,2) ,Xs (:,3) ) ; 
Xs (:, 1) ,mod(Xs (:,2)-2,L)+l,Xs (: ,3) ) ; 
Xs (:, 1) ,Xs{ :,2) ,mod (Xs { : , 3) -2, L) +1) ; 



zeros (12, lops, la) ; 
sros (length (Xc) , 12, lops) ; 
zeros (length (Xc) , 12, lops) ; 
kplaq-zeros (Nsweeps, la) ; 

at up X for operators 
X=Xc (1 :L"2) ; 
for r=l:12 

eval ( [' X' num2str (r) ' -repmat (X, [1 lops] ) ; ' ] ) 
nd % for r 

for k-l: lops 
c=ops (k) ; 



s (k, : ) )=ovrlaxJ(s 



end 

md % OR sweeps 



3val ( ['; 
%quit 



XI ( : 


, k) -xp (XI ( 


, k) , 1) ; 


X2 ( : 


, k) -Xp (Xz ( 


, k) , 1) ; 


X4 ( : 


, k) =Xm(X4 ( 


,k) , 1) ; 


X5 [ : 


, k) -Xm(X5 ( 


,k) , 1) ; 


X7{: 


,k) -Xm(X7 ( 


,k) ,2) ; 


X9{: 


,k)-Xp(X9( 


,k),2); 


XIO ( 


:,k)=Xp{X10 (:,k) ,2) ; 


X12 1 


: ,k) -Xm(X12 ( : , k) , 2) ; 



for hmc-l:PN 
%while hmc 

% hmc- hmc +1 
%echo off 

% compund sweeps to updateJ gauge fields Psi 

% loop over parity [Stagger) 

% loop over direction 

% calculate stapleJ sum R for link U 

% update J U 

for mu-l:D 

for k-1 :Str)ggpr 



for 1-1 :r 






XI 1 


-k) 


-Xm(Xl ( 


,k) ,2) 


X2 ( 


,k) 


=Xp (X2 ( 


,k) ,2) 


X4 ( 


, k) 


=Xp (X4 ( 


,k) ,2) 


X5 ( 


, k) 


=Xm(X5 ( 


,k) ,2) 


X7 ( 


,k) 


-Xp (X7 ( 


,k) , 1) 


X9 ( 


-k) 


-Xp (X9 ( 


,k) ,1) 



end 

for 1-1: a 
X3(:,k)- 
X6{:,k)- 
X8{:,k)- 
Xll [ : ,k) 

pnd 



Xp(X3(:,k),2); 
Xm(X6 ( : , k) , 2) ; 
Xp(X8(:,k),l); 
-Xm(Xll ( : , k) , 1) ; 



Figure C.3: Programs used for the lattice calculation. 
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for h=l : skip :N sweeps 
for a-4:la 

for pk-l:4 

f lag=0; 

eval ( [ ' load ' 



d3-real (K (sub2ind {sM, 2*X5 ( : , 
d4-imag(K (sub2ind {sM, 2*X5 ( : , 



infile num2str(h.j 



% calculate 0++ 

for t-l:L 

x-[l:L-2]'; 

xl-Xp (X, 1) ; % X + mu 

x2-Xp (X, 2J ; % X + nu 

x-(t-lt*L'-2+[l:L-2]' 

switch pk 

% no change 
disp (h) 



% change side 2 

Psi (K~2+ [1:N''2] ,x (1:LJ ) = 

% change both 

Psi (1 :N"2,x (1:L;L''2) ) = 

Psi (N-'2+ [1:N''2] ,x (1:LJ ) = 



Psi {1 :N'2,x (1 :L:L~2) ) ; 
* Psi (N''2+ [1 :N"2] , x{l :L) ) ; 



Psi {1 :N"2,x (1 :L;L''2) ) ; 
* Psi (N''2+ [1 :N"2] , x{l :L) ) ; 



for k-l:L''2 

M(2*k-l:2*k,2*xl (k)-l:2*xl (k) ) — alph{a) . 
M (2*xl (k) -1 : 2*xl (k) , 2*k-l : 2*k) =-alph (a) . 
M(2*k-l:2*k,2*x2 (kj -l;2*x2 (k) ) — alph{a) , 

M(2*x2{k)-l:2*x2(kJ,2*k-l:2*k)- 

end % for k 

M-spdiags{ones (sM (1) , 1) , 0, M) ; 
K-eye (sM) ; 

%K-inv[speye [size (M) ) -alph . *M) ; 
if isempty (Pc) 
Pc^colmmd (M) ; 

Pc-spar3e{l:sM(l) ,Pc,l,sM(l) ,sM(l) ) ; 
K^mexmysolve [M, K, Pc) ; 



'su2 (Psi{l:N'2,x(k) ) ) ; 
'su2 (Psi(l:N*2,x(k) ) ) ' ; 
'su2 (Psi (N''2+ [1 :N''2] ,x( 



-alphfa) .*su2 (Psi {N''2+ [1 :N^2 ] , 



bl.*c3 
l.*b3. 
b4.*c2 
b3.*cl 
b2,*c2 
b3.*c4 
b2,*c3 
b3.*cl 
4.*b4. 
S-'bl. 
4.*bl. 



0((t-l) 
-2,*al. 
*dl-a3.*bl.*cl.-d4 
-al.*b3 
^cl.*d3+2.*a2.-'b2.*. 

-2. *a4 . 
'■b2. -c2 
-2. -al 
'-b3. -c4 
+2. *a2 
'bl . 'c2 



*d4+2.*a4, 
*d2+2. 
*dl-2.*a2, 
*dl-2. 
*d4+2.*a2, 
*d4-2.*a3, 
c3 . *d4+al . 
c4 . *d3+a3 . 
c2 . *d3-a2 . 



-2. *a2. 
*b2.*c4. 

+2.*a4. 
*bl.*c3. 

-2.*al. 
*b2.*cl. 

-a2 . *bl 
*b4 . *c2 . 

+ a4 . ■'bl 
'b4 . "^4 . 

+ a4 . ■'b4 
■b4 . *cl 



2+k,2, s)= ... 

*b3.*c3.*dl+2.-al.*bl.*cl.*dl-2,*a3.*b3.*cl.*dl-2.'-a 
-al.*bl.*c3.-d4+a3.*b3.'c3.-d4 , . . 
. *cl . *d4-a3 . *bl . *cl . *d3-al . *bl . *c3 . *d3+a3 . *b3 . *c3 . 
c2 . *d2 + 2 . *a4 . ■*b4 . *c2 . *d2 + 2 . *a2 . ■'b4 . *c4 . *d2 . . . 
*b2 . ■'C4 . *d2+2 . *a2 . *b2 . -04 . 'd4+2 . 'a4 . *b4 . •c4 . 'd4-2 . '"a. 
*d4-2 . *al . 'bl . 'c2 . *d2+2 . *a3 . -bS. 'c2 . *d2 . . . 
*b2 . 'cl . *d2-2 . *a2 . *bl . -cl . 'd2+2 . 'a3 . *b4 . -cl . 'd2-2 . *a. 
*d2+2 . *a3 . 'bl . *c4 . *d2-2 . 'al . •b4 . 'c3 . 'd2 . . . 
*b3 . ■'c3. *d2-2 . *a3 . *b2 . -cS. 'd2-2 . 'a4 . *bl . -cS. 'd2-2 . 'al . 
*dl+2 . *a3 . ■'b4 . *c2 . *dl-2 . *a4 . ■'b3. 'c2 . *dl . . . 

a4 . *b4 . ■'cl . *dl-2 . *al . *b4 . *c4 . *dl+2 . *a2 . 
c4.*dl-2.*a2.*b4.*c3.*dl ... 
al.*bl.*c4.'d4+2.'a3.-b3.*c4.'d4+2.'-a 
c3.-d4+2.'a4.*b3.*c3.-d4 ... 
a3 . *bl . *c2 . *d4-2 . *al . *b4 . *cl . *d4+2 . '"a. 
cl . *d4-al . *b2 . *c4 . *d4 , . . 
. *c4 . *d4+a3 . *b4 . *c4 . *d4-a4 . *b3 . *c4 . *d4+a2 . *b2 . *c3 . '■^4 +a 
d4-a2 . ■*b3. *c2 . *d4+a3 . *b2 . *c2 . *d4 , . . 

*c2 . *d4-a2 . *b4 . "'cl . *d4+a4 . *b2 . *cl, *d4-al . *b2 . *c4 . '■d3-a 
■'d3-a4. 'b3 . *c4 . *d3+a2 . *b2 . *c3 . *d3 . . . 

. 'c^ . 'd3 + al . *b4 . ■'c2 . *d3-a2 . *b3 . *c2 . *d3+a3 . *b2 . *c2 . '"ai + a 
■'d3 + a4 . •b2 . *cl . *d3; 



■b2 . *cl .*dl-2 . 
■dl-2.*a4.*bl. 
■b2.*c3.*dl-2. 
■d4-2.'a3.*b4, 
■b3.*c2.*d4-2. 
■d4-2.*a4.*bl. 



al=real (K ( sub2ind (sM, 2'X-l , 2'X5 ( : , s j -1) ) ) ; 
a2=imag (K ( sub2ind [sM, 2*X-1 , 2*X5 ( : , s) -1) ) 1 ; 
a3-real (K ( sub2ind (sM, 2*X, 2*X5 ( : , s) -1) ) ) ; 
a4-imag(K(sub2ind(sM, 2*X,2*X5( :,s)-l) ) } ; 



bl-real (K (sub2ind (sM, 2*X5 ( : , 
b2-iiiiag (K (sub2ind (sM, 2*X5 ( : , 
b3-real (K [ sub2ind {sM, 2*X5 ( ; , 
b4-iiiiag (K ( sub2ind {sM, 2*X5 ( : , 

cl-real (K ( s-.;b2ind ( sM, 2-'X6 ( : , 
c2-iiiiag (K { s-.-,r..2irr.i ( :-.M, 2-'X6 ( : , 
c3-real (K ( s-^b2:.nci ( sM, 2-'X6 ( : , 
c4-iinag (K ( sub2ind (sM, 2*X6 ( : , 

dl-real (K (sub2ind {sM, 2'Xl ( : , 
d2-imag (K ( sub2ind [sM, 2*X1 ( : , 
d3=i:eal (K ( 3ub2ind ( sM, 2*X1 ( ; , 
d4=imag (K { s-jb2ind ( sM, 2*X1 ( ; , 



s)-l,2'X6(:,s)-l))); 
st-l,2*X6(;,s)-in); 
st,2*X6{:,s)-l))t; 
s) ,2*X6{ :,s)-l) ) ) ; 

3) -1,2*X1(: , sj -!))); 
3) -1,2*X1(: , 3) -1) ) ) ; 
3) ,2-Xl ( s)-l) ) j ; 
s) ,2*X1 ( :,s)-l) ) ) ; 

st-l,2'X-in); 
s)-l,2'X-in); 
s) ,2*X-1) ) ) ; 
s) ,2*X-1) ) ) ; 



O ( (t-1) *L-2+k, 3,3)= ... 

-2 . *al . *b3 . *c3. *dl+2 . *al . *bl . 'cl . *dl-2 
■cl . *d4-al . -bl . 'c3 . *d4 + a3 . ■•b3. 'c3 . *d4 
-al . *b3 . 'cl . *d4-a3 . *bl . -cl . 'd3-al . 'bl . ^ 
cl . *d3 + 2 . *a2 . ■'b2 . 'c2 . *d2 + 2 . *a4 . ■'b4 . 'c2 . *d2 + 2 . *a2 . 'l 
-2 . •a4 . *b2 . ■'C4 . *d2+2 . *a2 . *b2 . ■'c4 . *d4+2 . 
*d4+2 . *a4 . ■'b2 . *c2 . *d4-2 . *al . 'bl . *c2 . *d2+2 . *a3 . ■'b3. 

-2 . *al . *b2 . ■'cl . *d2-2 . *a2 . *bl . ■'cl . *d2+2 . 

■d2+2 . *a3 . *bl . *c4 . *d2-2 . *al . *b4 . 
■b3 . *c3 . *d2-2 , *a3 . *b2 . *c3 . *d2-2 , 
'd1 +2 . *a 3 . *b4 . *c2 . *dl-2 . *a4 . *b3 . 



|bl . *c3 . *dl-a3 . ■'bl . 

*b3. 
|b4.-*c2 
b3.*cl 
b2,*c2 
b3.*c4 
b2.*c3.*d4+2 



•a3 . *b3 . ■'cl . *. 



c3. 'd3+a3 . ■*! 



'a3.*b4.^'cl 



*d2+2.^ 



-dl-2.*a2.*bl . 
*dl-2 



';i2.*b2. 



*dl- 



*dl- 



dl-2 . 
b3 . 'c^ . 
'd4-2 . 
*d2-2 . 
*d2-2 . 
*dl+2 . 



al-real (K [sub2ind (sM, 2*X-1, 2*X1 ( : , s) -1) ) ) ; 
a2-imag (K [sub2ind(sM, 2*X-1,2*X1 ( : , s) -1) ) ) ; 
a3-real (K (3ub2ind (sM, 2*X, 2*X1 ( ; , s) -It) ) ; 
a4-imag (K (sub2ind (sM, 2*X, 2*X1 ( : , s) -1) ) ) ; 



bl-real (K (3ub2ind (sM, 2*X1 ( : , 
b2-imag (K (3iib2ind (sM, 2*X1 ( : 
b3-real (K {3ub2ind (sM, 2*X1 ( : 
b4-imag (K {sub2ind (sM, 2*X1 ( : 

cl=real (K (sub2ind (sM, 2'X2 ( : 
c2 = imag {K (subzi.Td (sM, 2*X2 ( ; 
c3=real (K (sub2ind (sM, 2*X2 ( : 
c4-iraag (K [sub2ind(sM, 2*X2 ( : 

dl-real (K {3iib2ind (sM, 2*X3 ( : , 
d2-imag (K (sub2ind (sM, 2*X3 ( : , 
d3-real (K (sub2ind (sM, 2*X3 ( : 
d4=imag (K (sub2ind (sM, 2*X3 ( : 



s)-l,2*X2(:,s)-l))); 
s)-l,2*X2 (:,s)-l) ) ) ; 
s) ,2*X2 (;, st-lt ) ) ; 
s) ,2*X2 (:, s) -1) ) ) ; 

3) -1, 2-X3 (: , 3) -1) j ) ; 
3) -1, 2*X3 ( : , 3) -1) j ) ; 
3) ,2*X3 (: , 3) -1) ) ) ; 
3) ,2*X3 (:, s) -1) ) ) ; 

s)-l,2*X-l)n; 
s)-l,2*X-l)n; 
s) ,2*X-1) ) ) ; 
s) ,2*X-lt ) ) ; 



ol . ■'c3 

1. *b3.^ 
□4.^'c2. 
o3.*cl. 

02. *c2. 

03. *c4. 

02. *c3. 

03. *cl. 

4. *b4.' 
2 . *bl . ■ 
4.*bl.' 



0( 

*dl-a3 . *bl 
cl. *d3-^2. '. 
*d4+2 . *a4 
*d2+2.*al 
*dl-2,*a2 
*dl-2,*a3 
*d4+2.*a2 
*d4-2.*a3 
c3.*d4+al 
c4.*d3+a3 
c2,*d3-a2 



(t-1 j *L"2+k, 1, s) = ... 
-2 . - al . *b3 . *c3 . •dl + 2 . ■'al . *bl . *cl . *( 
. ■'cl . 'd4-al . *bl . ■'c3. 'd4+a3 . *b3 . ■'c3. 

-al . *b3 . *cl . ■'d4-a3. 'bl . 'cl . ■'d3-al . 
a2 . 'b2 . ■'c2 . 'd2+2 . 'a4 . 'b4 . ■'c2 . 'd2+2 . 

. *b2 . *c4 . ■'d2+2 . ■'a2 . *b2 . *c4 . *( 



b2 . ■'c2 . 'd4-2 . *al . *bl . 



■'al. *b2.^ 



*d2-2 . 



b3.*c4.*d2+2.*a3. 
+2.*a2,*b3.-c3.*d2-2. 
bl.*c2,*dl+2,*a3.*b 
-2.*a2.-*b2.*cl.*dl- 
b2 . ■'C4 . 'dl-2 . 'a4 . 'fc 
+ 2 . 'a4 . 'b2 . 'c3 . ■dl-2 . 
bl . •c3 . 'd4-2 . 'a3 . 



^c2.*d2+2.*a3.*b3 
^a2.*bl.*cl.*d2+2 
'c4.*d2-2.*al.*b4 
'a3.*b2.-c3.*d2-2 
'c2,*dl-2,*a4.*b3 
^a4.*b4.'cl.^'dl-2 
■c4 . 'dl-2 . 'a2 . 'b4 



-2 . -al . 'b3 . 'c2 . 
b2 . ■'cl . 'd4-2 . 'a. 



4-2. 



. 'cl . 'd4-al . 'b2 . '< 
-a2 . *bl . 'c4 . ■d4+a3 . 'b4 . 'c4 . 'd4-a4 . 'b3 
b4 . ■'c2 . 'd4-a2 . 'b3 . ■c2 . 'd4+a3 . 'b2 . •c2 . '■ 
+a4 . *bl . 'c2 . *d4-a2 . 'b4 . 'cl . ■'d4+a4 . *b2 . 
b4 . *c4 . *d3-a4 . *b3 . *c4 . *d3+a2 . *b2 . *c3 . *. 
+a4 . *b4 . *c3 . *d3+al . *b4 . *c2 . *d3-a2 , *b3 . 
b4 . *cl . *d3+a4 . *b2 . *cl . *d3; 

al-real (K(3iib2ind(sM, 2'X-l, 2*X3 ( : , s) -1) ) ) ; 
a2-imag (K (3iib2ind (sM, 2'X-l, 2*X3 ( : , s) -1) ) ) ; 
a3-real (K (sub2ind (sM, 2'X, 2'X3 ( : , si -1 1 ) ) ; 
a4=imag (K (sub2ind(sM, 2'X, 2'X3 ( : , si -1) ) ) ; 



.*a3.*b3.*cl.*dl-2. 
'c3.^'d3+a3.*b3.*c3. 

*d4-2. 
*d2-2. 
*d2-2 . 
*dl+2. 
■d4+2 . 
'd4+2 . 

'c4 . ■d4+a2 . 'b2 . 'c3 . 
4 ... 

cl . *d4-al . 'b2 . 'c4 . 
3 . . . 

c2.*d3+a3.*b2.*c2. 



'a4. *b4.^ 
'c2.*d2 
'a3. *b4.^ 
►c3.*d2 
'a4,*bl.^ 
'c2,*dl 
'al. *b4 
. •c3. 'dl 
'a3 . 'b3.^ 
'c3 . 'd4 



bl=real (K (sub2ind (sM, 2'X3 ( : 
b2=imag (K (sub2ind(sM, 2'X3 ( : 
b3=real (K (sub2ind (sM, 2'X3 ( : 
b4 = imag (K (sub2ind(sM, 2'X3 ( : 

cl-real (K (sub2ind (sM, 2*X4 ( : 
c2-imag (K {3ub2ind (sM, 2*X4 ( : , 
c3-real (K (3ub2ind (sM, 2*X4 ( : , 
c4-imag (K (3ub2ind (sM, 2*X4 ( : , 

dl-real (K {3ub2ind (sM, 2*X5 ( : 
d2^-imag (K(Biih2i nd (sM, ^'X.ii ( : 



3l -1, 2'X4 ( : , 3l -11 1 ) ; 
3l -1, 2'X4 ( : , 3l -11 1 ) ; 
3l ,2'X4 (: , 3l -11)1; 
3l ,2'X4 (: , 3l -11)1; 

3l-l,2*X5 (: ,sl-ll 1 ) ; 
s)-l,2*X5(:,s)-l)l); 
s),2'X5(:,s)-in); 
s),2'X5(;,st-lt)); 

s)-l,2*X-l) ) ) ; 



*bl . *c2 . 'd3-a2 



■b4 . 'c4 . 'd3-a4 . 'b3 . 'c4 . 'd3+a2 . *b2 . * 
+a4 . 'b4 . 'ci . ^d3+al . 'b4 . *c2 . *d3-a2 . 
■b4 . 'cl . 'd3+a4 . *b2 . 'cl . 'd3; 



■b3 . *c2 , *d3+a3 . *b2 . *c2 . 'a 



al-real (K ( sub2ind (sM, 2*X-1 , 2*X6 ( : , s) -1) ) ) ; 
a2-ii[iag(K(sub2ind[sM,2*X-l,2*X6(:,sJ-lt ) ) ; 
a3-real (K ( sub2ind (sM, 2*X, 2*X6 ( : , s) -1) ) ) ; 
a4-imag(K(sub2ind(3M, 2*X,2*X6( :,s)-l) ) ) ; 



bl-real (K (sub2ind (sM, 2*X6 ( : , 
b2-imag (K (sub2ind (sM, 2*X6 ( ; , 
b3-real (K (sub2ind (sM, 2*X6 ( ; , 
b4-iiiiag (K (sub2ind (sM, 2*X6 ( : , 

cl-real (K ( sub2ind [sM, 2*X1 ( : , 
c2-ii!iag (K ( sub2ind [sM, 2'Xl ( : , 
c3-real (K ( sub2ind [sM, 2*X1 ( : , 
c4-iinag (K (sub2ind {sM, 2*X1 ( : , 

dl-real (K (sub2ind (sM, 2'X2 ( : , 
d2-imag (K ( sub2ind [sM, 2*X2 ( : , 
d3=i:eal (K ( 3ub2ind (sM, 2'X2 ( ; , 
d4=imag (K { s^jb2ind (sM, 2'X2 ( ; , 



s)-l,2'Xl(:,s)-l))); 
s)-l,2*Xl(;,s)-lt)); 
St ,2*X1{ :,s)-l) ) ) ; 
si ,2*X1 { :,s)-l) 1 1 ; 

3l -1,2'X2 (: , 3l -11)1; 
3l -1,2'X2 (: , 3l -11 ) 1 ; 
si ,2'X2 ( s)-l) 1 1 ; 
si ,2*X2 ( :,s)-l) 1 1 ; 

st-l,2'X-in); 
s)-l,2"X-lt ) ) ; 
s),2*X-l))t; 
s) ,2*X-1) ) ) ; 



bl.-*c3. 

1. *b3. ' 
b4.-*c2. 
b3.*cl. 
b2,*c2, 

03. *c4, 

02. *c3. 

33. *cl. 

4 . ■'b4 .' 
2.*bl.■ 
4.*bl.■ 



*dl-a3 . 'bl 
cl . 'd3 + 2 . ', 
'd4+2 . 'a4 . 
*d2+2.*al. 
-dl-2.*a2. 
*dl-2.'a3. 
*d4+2.*a2. 
*d4-2.*a3. 
c3.*d4+al. 
c4.*d3+a3. 
c2.*d3-a2. 



O ( (t-1] 'L-2 + k, 4, 3l = ... 

-2 . *al . 'b3 . 'c3 . 'dl+2 . 'al . 'bl . *cl . 'dl-2 
1 . -bl . 'c3 . 'd4+a3 . •b3. 'c3 . 'd4 
•cl . *d4-a3 . 'bl . 'cl . 'd3-al . 'bl . •c3. '■ 
c2 . 'd2+2 . 'a4 . *b4 . 'c2 . 'd2+2 . 'a2 . *b4 . '■ 



*cl . 'd4-, 
-al . 'b3 . 
.2.'b2. 'c: 
-2. *a4 . 
b2. 'c2. 
-2. *al. 
b3. *c4. 



*d2-2. 'a2 



'b2 . *c4 . 'd4+2 . 'a4 
'd2+2 . 'a3 . *b3. 'c2 
'bl.^'cl.'d2+2.'a3 



'dl-2 . 



'd4+2 . 



'd4+2 . 



b2 . 'cl . 
-a2 . 'bl 
b4 . 'c2 
■fa4 .'bl 



b4.*cl 



d2+2.*a3.*bl. *c4 
b3.*c3.*d2-2,*a3 
■dl+2.*a3.*b4,*c2 
b2.^'cl.*dl-2. *a4 
dl-2 . 'a4 . -bl . 'c4 

d4-2 . 'a3 . 'b4 . 'c3 

b3 . 'c2 . ^d4-2 . 'a3 

d4-2 . 'a4 . 'bl . 'cl 

'C4.'d4+a3.'b4.'c4.'d4-a4.'b3.'c4 
'd4-a2 . ■b3 . 'c2 . 'd4+a3 . •b2 . 'c2 . 'd4 

*c2 . *d4-a2 . 'b4 . 'cl . 'd4+a4 . 'b2 . 'cl 
*d3-a4 . *b3 . *c4 . *d3+a2 . *b2 . *c3 . *d3 
. *c3 . *d3+al . *b4 . *c2 . *d3-a2 . *b3 . *c2 
*d3+a4.*b2.*cl.*d3; 



b3 . *cl . 'dl-2 . 

d3+a3 . 'b3 . -^3 . 
d2 ... 
*c4 . 'd4-2 . 

'cl. 'd2-2 . 

*c3.*d2-2 . 

*c4 . *dl+2 . 
'dl . . . 

'b3. 'c4 . '■■.-■.4(2. 
'd4 . . . 

. 'a4+2 . 
d4+a2 . 'b2 . 'c3 . 
d4-al.'b2.'c4. 
d3+a3.*b2.*c2 . 



al-real (K (sub2ind (sM, 2*X-1 , 2*X2 ( : , s) -1) ) ) ; 
a2-iirag(K (sub2ind(sM, 2*X-1, 2*X2 ( ; , s) -11 ) 1 ; 
a3-real (K ( sub2ind (sM, 2'X, 2-'X2 ( : , s) -1) 1 1 ; 

a4=imag (K ( 3ub2ind (sM, 2'X, 2*X2 ( : , s) -1) 1 ) ; 



bl = real (K { s^jb2ind (sM, 2'X2 ( : , 
b2=imag (K { 3ub2ind (sM, 2'X2 ( : , 
b3=real (K ( 3ub2ind (sM, 2'X2 ( : , 
b4=imag (K ( 3iib2ind (sM, 2'X2 ( : , 

cl-real (K ( sub2ind (sM, 2*X3 ( : , 
c2-imag (K (sub2ind (sM, 2*X3 ( : , 
c3-real (K (sub2ind (sM, 2-X3 ( : , 
c4-imag (K (sub2ind (sM, 2*X3 ( : , 

dl-real (K (sub2ind (sM, 2*X4 ( ; , 
d2^i mag (K ( snh2 -j nd ( ■■gM, ? 'X4 ( : 



3l -1, 2'X3 ( : , s) -1) ) 1 ; 
3l -1, 2'X3 ( : , 3l -11 ) 1 ; 
3l ,2'X3 ( s) -1) 1 1 ; 
3l ,2'X3 ( s)-l) 1 1 ; 

si -1,2*X4 (:, si -11)1; 
si -1,2*X4 (:,sl -11 ) ) ; 
s) ,2*X4 ( :,s)-l) ) 1 ; 
St ,2*X4 ( :,s)-l) ) 1 ; 

s)-l,2*X-lt ) ) ; 



Figure C.4: Programs used for the lattice calculation. 
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ol.*c3 

1. *b3. 
o4.*c2 
o3.*cl 

02. *c2 

03. *c4 
32.*c3 

03 . -cl 

4 . *b4 . 
2.'bl. 
4 . *bl . 



d3-real (K {3ub2ind (sM, 2*X4 ( : , s) ,2*X-1) ) ) ; 
d4-imag (K {3iib2ind (sM, 2*X4 ( : , s) , 2*X-1) ) ) ; 

0( (t-l)*L'-2+k,5,3)= . . . 

-2 . *al . *b3. -c3 . *dl+2 . *al . *bl . *cl . *dl-2 . *a3 . *b3 . -cl . *dl-2 . 
. *dl-a3 . -bl . *cl . *d4-al . *bl . *c3 . *d4+a3 . -b3 . *c3 . *d4 ... 

-al .*b3 . *cl . *d4-a3. *bl . *cl . -dS-al . *bl . *c3 . *d3+a3. *b3 . *c3 . 
'cl,*d3+2 . *a2 . *b2 . *c2 . *d2+2 . 'a4 . *b4 . ■'c2 . *d2+2 

-2 . ■*a4 . 'b2 . *c4 . *d2 + 2 
.*d4+2 .*a4 . *b2 . ■*c2 . *d4-2 . *al . ■'bl 

-2 . ■'al . *b2 . *cl . ■d2-2 
.*d2+2 .*al . *b3 . *c1 . 'd2+2 . 'a3 . *bl 

+2.^'a2.*b3.*c3.^*d2-2 
.*dl-2 .*a2 . *bl . *c2 . 'dl+2 . 'a3 . *b4 



*a2 . *b2 . ^ 



*dl-2 



'dl-2. 



*a2 . *b2 . ■'C4 . ■'d4 + 2 . ■'a4 . 
"cl . *d2+2 . *a3 . ■'b3 .'c2. 
■a2 . *bl . *cl . ■'d2 + 2 . ■'a3 . 
■'C4 . 'd2-2 . 'al . *b4 . ■'c3. 
■*a3.-b2.*c3.^*d2-2.^*a4. 
■'c2 . 'dl-2 . 'a4 . *b3 . ■'c2 . 
*a4 . *b4 . *cl . *dl-2 . *al . 
*c4.*dl-2.*a2.*b4.*c3. 
*al.*bl.-c4.*d4+2.*a3. 
*c3.*d4+2,*a4.*b3.*c3. 
*a3.*bl.*c2.*d4-2.*al. 
*cl. *d4-al . *b2 . *c4 , *d4 



b2.*c4.*dl-2.*a4.*bl 
+2 . *a4 . *b2 . -c3 . *dl-2 
*d4+2 . *a2 . *bl . *c3 . *d4-2 . *a3 . *b4 
-2.*al.*b3.*c2.*d4-2 
*d4-2 . *a3 . *b2 . ■*cl. *d4-2 . *a4 . *bl 

-a2 . 'bl . *c4 . *d4+a3. *b4 . *c4 . *d4-a4 , *b3 . *c4 
c3. *d4+al . *b4 . ■'c2 . *d4-a2 . *b3 . *c2 . *d4+a3 . *b2 . *c2 , *d4 

+ a4 . *bl . *c2 . ■'d4-a2 . *b4 .'cl. ■'d4 + a4 , *b2 . *cl 
c4 . *d3 + a3 . *b4 . ■•c4 . 'd3-a4 . *b3 .'ci. 'd3+a2 . *b2 . *c3. *d3 

+ a4 . *b4 . 'c3 . ■'d3 + al . 'b4 . *c2 . ■'d3-a2 . *b3 . *c2 
c2 . *d3-a2 . *b4 . ■'cl . 'd3+a4 . *b2 . ■'cl . 'd3; 

al-real (K (sub2ind (sM, 2'X-l, 2'X4 ( : , s) -1) ) ) ; 
a2-imag (K [sub2±nd(sM, 2*X-1, 2*X4 ( : , s) -1) ) ) ; 
a3-real (K (sub2ind (sM, 2*X, 2*X4 ( : , s) -1) ) ) ; 
a4-imag (K (sub2ind (sM, 2*X, 2*X4 ( : , s) -1) ) ) ; 



'd2-2. 
•d2-2. 
*dl+2 . 
*d4+2. 
*d4+2. 



^dl 
^b3. 
^d4 

^b4.'cl 
d4+a2,*b2.*c3 
d4-al.*b2.*c4 
d3+a3.*b2.*c2 



bl-real (K (3ub2ind (sM, 2*X4 ( 
b2-imag (K {3iib2ind (sM, 2*X4 ( 
b3-real (K {3iib2ind (sM, 2*X4 ( 
b4-imag (K {sub2ind (sM, 2*X4 ( 

cl=real (K (sub2ind (sM, 2'X5 ( 
c2 = imag {K (subzi.Td (sM, 2*X5 ( 
c3=real (K (sub2ind (sM, 2'X5 ( 
c4-iraag (K [sub2±nd(sM, 2*X5 ( 



s)-l,2'X5(:,s)-l)t); 
s)-l,2*X5 (:,s)-l) ) ) ; 
s) ,2*X5(;, s)-lt ) ) ; 
s) ,2*X5(:, s) -1) ) ) ; 

3) -1, 2-X6 (: , 3) -1) j ) ; 
3) -1, 2'X6 ( : , 3) -1) j ) ; 
3) ,2'X6(: , sj -1) ) 1 ; 
3) ,2*X6(:, s) -1) ) ) ; 



dl-real (K {3ub2ind (sM, 2*X6 ( : , s) -1, 2*X-1) ) ) ; 



d2-imag (K (sub2ind (sM, 2*X6 ( 
d3-real (K (sub2±nd (sM, 2*X6 ( 
d4=imag (K (sub2ind (sM, 2*X6 ( 



s)-i,2*x-i)t); 
•x-it)) ; 
*x-it ) ) ; 



( (t-1 j *L"2+k, 6, s) = ... 

-2 . ■'al . *b3 . 'c3 . ■'dl + 2 . ■'al . *bl . *cl . 
'c3. *dl-a3 . *bl . *cl . *d4-al . *bl . ■'c3. 'd4+a3 . *b3 . ■'c3 

-al . *b3 . 'cl . ■'d4-a3. 'bl . 'cl . •d3-al . *bl . *c3 . 
33 . ■'cl . ■*d3+2 . *a2 . *b2 . ■'c2 . *d2+2 . *a4 . 'b4 . ■'c2 . 'd2+2 . *a2 . *b4 . 

-2 . ■'a4 . *b2 . *c4 . *d2+2 . *a2 . *b2 . *c4 . *d4+2 . ■'a4 
'c2 . *d4+2 . *a4 . *b2 . *c2 . *d4-2 . *al . *bl . *c2 . *d2+2 . *a3 . *b3 . *c2 

-2 . *al . *b2 . *cl . *d2-2 . *a2 . *bl . *cl . *d2+2 . *a3 
'cl . *d2+2 . *al . *b3 . *c4 . *d2+2 . *a3 . *bl . *c4 . *d2-2 . *al . *b4 . *c3 

+2 . *a2 . *b3 . *-c3 . *d2-2 . *a3 . *b2 . *c3 . *d2-2 . *a4 
'c2.*d1-?.*a2.*b1 . *c? . *d1 +2 . *a3 . *b4 . *c? . *dl -2 . *a4 . *b3 . *c2 



*b2.*cl.*dl- 
'(II 



dl-2.*a3.*b3.*cl.*dl-2.*a3. 



d3+a3.*b3.*c3.*d3- 
c4.*d2 ... 
*b4.*c4.*d4-2.*a2. 
*d2 ... 

*b4.*cl.*d2-2.*a4. 
*d2 ... 

*bl.*c3.*d2-2.*al. 
*d1 ... 

*b4.'c4.^'dl + 2.^'a2. 



bl.*c3 

1. ^'b3.^ 
34 . *c2 

03. *cl 
d2. *c2 
b3.*c4 
b2,*c3 
b3.*cl 

4. *b4. 

2. ^'bl. 
4.*bl. 



bl . *c3 . 'dl-a3 . ■'bl . 

■'b3. 
b4.-*c2 
b3.*cl 
b2,*c2 



d3-real (K (sub2ind {sM, 2*X10 { : 
d4-imag(K (sub2ind {sM, 2*X10 { : 



0((t-l) 

-2,*al. 
-dl-a3.*bl.*cl.-d4 
-al.*b3 
^cl.-d3+2.*a2.^'b2.*. 



2+k,; 



*d4+2.*a4 
*d2+2.*al 
*dl-2.*a2 
'dl-2.*a3 
*d4+2.*a2 
*d4-2.*a3 
c3 . 'd4+al 
c4 . 'd3+a3 
c2 . 'd3-a2 



-2. *a4 
■'b2. 'c2 

-2. -al 
■'b3. 'c4 

+2. *a2 
■'bl . 'c2 



-2. *a2. 
*b2.*c4. 

+2.*a4. 
*bl.*c3. 

-2.*al. 
*b2.*cl. 

-a2 . *bl 
'■b4 . 'c2 . 

+a4 . 'bl 
'b4 . 'c4 . 

+a4 . 'b4 
■b4 . 'cl 



*b3.*c3.*dl+2,*al.*bl.*cl.*dl-2,*a3.*b3.*cl.*dl-2.'a3 
-al . *bl . *c3 . -d4+a3 . *b3 . 'c3 . -d4 ... 

. *cl . ■'d4-a3 . *bl . *cl . *d3-al . *bl . *c3 . *d3+a3 . *b3 . *c3 . 'd3 
c2 . 'd2 + 2 . 'a4 . ■'b4 . *c2 . 'd2 + 2 . 'a2 . ■'b4 . 'c4 . 'd2 . . . 
'b2 . ■'C4 . *d2+2 . 'a2 . 'b2 . •c4 . 'd4+2 . 'a4 . 'b4 . •c4 . 'd4-2 . 'a2 
'd4-2 . 'al . *bl . 'c2 . 'd2+2 . 'a3 . 'bS. 'c2 . 'd2 . . . 
'b2 . *cl . *d2-2 . 'a2 . 'bl . 'cl . 'd2+2 . 'a3 . 'b4 . 'cl . 'd2-2 . 'a4 
'd2+2 . 'a3 . *bl . 'c4 . 'd2-2 . 'al . •b4 . 'c3 . 'd2 . . . 
'b3 . *c3. *d2-2 . 'a3 . 'b2 . 'cS. 'd2-2 . 'a4 . 'bl . 'cS. 'd2-2 . 'al 
'dl+2 . 'a3 . *b4 . *c2 . 'dl-2 . 'a4 . *b3. 'c2 . 'dl . . . 
*b2 . *cl .*dl-2 . *a4 . *b4 . ■'cl . *dl-2 . *al . *b4 . *c4 . *dl+2 . 'a2 
*dl-2.*a4.*bl.*c4.*dl-2.*a2.*b4.*c3.*dl . . . 
*b2.*c3.*dl-2.'al.*bl.*c4.'d4+2.'a3.-b3.*c4.'d4+2.'a: 
*d4-2.'a3.*b4,'c3.-d4+2.'a4.*b3.*c3.-d4 , . . 
*b3 . *c2 . *d4-2 . *a3 . *bl . *c2 . *d4-2 . *al . *b4 . *cl . *d4+2 . 'a2 
*d4-2 . *a4 . *bl . *cl . *d4-al . *b2 . *c4 . *d4 ... 
. *c4 . *d4+a3 . *b4 . *c4 . *d4-a4 . *b3 . *c4 . *d4+a2 . *b2 . *c3 . 'd4 + 
'd4-a2 . ■'b3. *c2 . *d4+a3 . *b2 . *c2 . *d4 , . . 

. *c2 . 'd4-a2 . 'b4 . ■'cl . *d4+a4 . *b2 . *cl, *d4-al . *b2 . *c4 . '63- 
'd3-a4. 'b3 . 'c4 . 'd3+a2 . *b2 . *c3 . *d3 . . . 

. ■c3 . 'd3+al . 'b4 . *c2 . *d3-a2 . *b3. *c2 . *d3+a3 . *b2 . *c2 . 'o3+ 
'd3+a4 . •b2 . 'cl . 'd3; 



al=real (K ( sub2ind (3M, 2'X-l , 2'XIO ( : , s) -1) ) j ; 
a2=imag (K ( sub2ind [sM, 2'X-l , 2'XIO [ : , s) -1) ) ) ; 
a3-real (K ( sub2ind (sM, 2*X, 2*X10 ( : , s ) -1 ) ) ) ; 
a4-imag(K(sub2ind(sM, 2*X,2*X10 (:,s)-l) } ) ; 



bl-real (K (sub2ind (sM, 2*X10 ( 
b2-iiiiag (K (sub2ind {sM, 2*X10 { 
b3-real (K [ sub2ind {sM, 2*X10 { 
b4-iiiiag (K ( sub2ind {sM, 2*X10 { 

cl-real (K ( sub2ind {sM, 2*X11 ( 
c2-iiiiag (K (sub2ind {sM, 2*X11 ( 
c3-real (K (sub2ind {sM, 2*X11 [ 
c4-iinag (K (sub2ind {sM, 2*X11 ( 

dl-real (K (sub2ind {sM, 2*X12 { 
d2-imag (K ( sub2ind [sM, 2*X12 ( 
d3=i:eal (K ( 3ub2ind ( 3M, 2'X12 ( 
d4=imag (K { 3-jb2ind ( sM, 2'X12 { 



i)-l,2*Xll (:,s)-l) ) ) ; 
i)-l,2*Xll (:,s)-l) ) ) ; 
i) ,2*X11 (;, s) -It) ) ; 
,2*X11 (:, s) -1) ) ) ; 

i)-l,2'X12 (: , 3) -1) j ) ; 
0-l,2'X12 ( : , 3) -1) j ) ; 
i) , 2'X12 (: , 3) -1) ) 1 ; 
i) ,2*X12 (:, s) -1) ) ) ; 

i)-l,2*X-l))t; 
i)-l,2*X-l) ) } ; 
i) ,2*X-1) ) ) ; 
i) ,2*X-1) ) ) ; 



O ( (t-1) 'L-2+k, 9,3)= ... 

-2 . *al . 'b3 . 'c3 . 'dl+2 . 'al . 'bl . *cl . 'dl-2 
■cl . 'd4-al . -bl . 'c3 . 'd4+a3 . -bS. 'c3 . 'd4 
-al . 'b3 . 'cl . *d4-a3 . 'bl . -cl . 'd3-al . 'bl 
cl . 'd3+2 . 'a2 . *b2 . 'c2 . 'd2+2 . 'a4 . *b4 . 'c2 . 'd2+2 . 'a2 . *1 
-2 . •a4 . 'b2 . *c4 . *d2+2 . 'a2 . 'b2 . *c4 . 'd4+2 . 
'd4+2 . 'a4 . *b2 . 'c2 . 'd4-2 . 'al . *bl . 'c2 . 'd2+2 . 'a3 . *b3. 

-2 . *al . *b2 . ■'cl . *d2-2 . *a2 . 'bl . ■'cl . 'd2+2 . 

■d2+2 . *a3 . *bl . *c4 . *d2-2 . *al . *b4 . 
■b3 . *c3 . *d2-2 , *a3 . *b2 . *c3 . *d2-2 , 
'dl +2 . *a 3 . *b4 . *c2 . *dl-2 . *a4 . *b3 . 



•a3 . 'b3 . *cl . '. 



•c3. 'd3+a3 . '1 



'a3.*b4.^'cl 



*d2+2.^ 



*d1-2.*a2.*bl . 



'c4.*dl-2.'a3. 
*c3.*d4+2.*a2. 



';i2.*b2. 



*dl- 



dl-2 
b3 . 'c^ 
'd4-2 
*d2-2 
*d2-2 
I 2 



*d3+a3 . 'b4 . •c4 . 'd3-a4 . 'b3 . *c4 . 'd3+a2 . 'b2 . *c3. *d3 . , . 

+a4 . *b4 . 'c3 . *d3+al . 'b4 . 'c2 . *d3-a2 . *b3 . *c2 . *d3+a3. *b2 . *c2 . *d3+a 
*d3-a2 . 'b4 . *cl . 'd3+a4 . 'b2 . *cl . 'd3; 

al-real (K (sub2±nd (sM, 2'X-l, 2*X12 ( : , s) -1) ) ) ; 
a2-imag (K [sub2±nd(sM, 2*X-1,2*X12 ( : , s) -1) ) ) ; 
a3-real (K (3ub2ind (sM, 2*X, 2*X12 { : , s) -1) ) ) ; 
a4-imag (K (3ub2ind (sM, 2*X, 2*X12 ( : , s) -1) ) ) ; 

bl-real (K (3ub2ind (sM, 2*X12 ( : , s) -1, 2*X7 ( : , s) -1) ) ) ; 
b2-imag (K (3iib2ind (sM, 2*X12 ( ; , s) -1, 2*X7 ( ; , s) -1) ) ) ; 
b3-real (K {3ub2ind (sM, 2*X12 ( ; , s) , 2*X7 { :,s)-l) ) ) ; 
b4-imag (K {sub2ind (sM, 2*X12 ( ; , s) , 2*X7 { : , s) -1) ) ) ; 



• u4 . 'd3 I u_i . 'bl . ' . ' Jj a4 . 'bJ . ' u4 . 'd3 I u2 . *b2 . * 
+a4 . 'b4 . 'ci . ^d3+al . 'b4 . *c2 . *d3-a2 . 
*c2 . 'd3-a2 . •b4 . 'cl . 'd3+a4 . *b2 . 'cl . 'd3; 



■b3 . *c2 , *d3+a3 . *b2 . *c2 . '0:3 + a 



al-real (K ( sub2ind (sM, 2*X-1 , 2*X7 ( : , s) -1) ) ) ; 
a2-ii[iag (K ( sub2ind (sM, 2*X-1 , 2*X7 ( : ,s) -1) ) ) ; 
a3-real (K ( sub2ind {3M, 2*X, 2*X7 { :,s)-l) ) ) ; 
a4-imag(K(sub2ind(3M, 2*X,2*X7 ( :,s)-l) ) ) ; 



cl=real (K (sub2ind (sM, 2'X7 ( 
c2 = imag {K (subzi.Td (sM, 2'X7 ( 
c3=real (K (sub2ind (sM, 2'X7 ( 
c4-iraag (K [sub2±nd(sM, 2'X7 ( 

dl-real (K {3iib2ind (sM, 2*X8 ( 
d2-imag (K (sub2ind (sM, 2*X8 ( 
d3-real (K (sub2±nd (sM, 2*X8 ( 
d4=imag (K (sub2ind (sM, 2*X8 ( 



3) -1) ) ) ; 



3) -1, 2'X8 ( : 
s) -1, 2'X3 ( ; 
3) , 2'X8 (: , 3) -1) ) ) ; 
3) ,2'XS (: , s) -1) ) 1 ; 

s)-l,2*X-l))); 
s)-l,2*X-l))); 
s),2*X-l))); 
s) ,2*X-1) ) ) ; 



bl-real (K (sub2ind (3M, 2*X7 ( : 
b2-imag (K (sub2ind {3M, 2*X7 ( ; 
b3-real (K (sub2ind {3M, 2*X7 ( ; 
b4-iiiiag (K (sub2ind {3M, 2*X7 ( ; 

cl-real (K ( sub2ind (sM, 2*X8 ( 
c2-ii!iag (K ( sub2ind (sM, 2*X8 ( 
c3-real (K ( sub2ind (sM, 2*X8 ( 
c4-iinag (K (sub2ind {3M, 2*X8 ( 

dl-real (K (sub2ind {3M, 2'X9 ( 
d2-imag (K ( sub2ind [sM, 2*X9 ( 
d3=i:eal (K ( sub2ind (sM, 2'X9 ( 
d4=imag (K ( 3ub2ind (sM, 2'X9 ( 



s)-l,2'X8(:,s)-l))); 
st-l,2*X8(;,s)-lt)); 
s) ,2*X8{ :,s)-l) ) ) ; 
s) ,2*X8{ :,s)-l) ) ) ; 

s) -1,2*X9(: , 3) -1) ) 1 ; 
s) -1,2*X9(: , 3) -1) ) 1 ; 
s) ,2*X9( :, s)-l) ) ) ; 
s) ,2*X9{ :,s)-l) ) ) ; 

st-1, 2'X-l))); 
s)-l, 2'X-l} ) ) ; 
s),2*X-l))}; 
s) ,2*X-1) ) ) ; 



ol . *c3 

1. 'b3.^ 
□4.^'c2. 
o3.*cl. 

02. *c2. 

03. *c4. 

02. *c3. 

03. *cl. 

4. *b4.' 
2 . *bl . ■ 
4.*bl.' 



0( 

*dl-a3 . 'bl 
cl. *d3+2. '. 
*d4+2 . *a4 
*d2+2.*al 
*dl-2,*a2 
*dl-2,*a3 
*d4+2.*a2 
*d4-2.*a3 
c3.*d4+al 
c4.*d3+a3 
c2,*d3-a2 



(t-1 j *L"2+k, 7, s) = ... 
-2 . - al . *b3 . 'c3 . •dl + 2 . ■'al . 'bl . 'cl . *( 
. *cl . 'd4-al . 'bl . *c3. 'd4+a3 . 'b3 . *c3. 

-al . *b3 . 'cl . *d4-a3. 'bl . 'cl . *d3-al . 
a2 . 'b2 . *c2 . 'd2+2 . *a4 . 'b4 . *c2 . *d2+2 . 

. *b2 . *c4 . ■'d2+2 . ■'a2 . *b2 . *c4 . *( 



b2 . ■'c2 . *d4-2 . *al . *bl 



■'al. *b2.^ 



*d2-2 



b3.*c4.*d2+2.*a3.*bl 
+2.*a2,*b3.-c3.*d2-2 
bl.*c2,*dl+2,*a3.*b4 
-2.*a2.-*b2.*cl.*dl-2 
b2 . *c4 . 'dl-2 . 'a4 . 'bl 
+2 . 'a4 . 'b2 . 'c3 . ■dl-2 
bl . •c3 . 'd4-2 . 'a3 . 'b4 



-2 . -al . 'b3 . 'c2 . 
b2 . *cl . 'd4-2 . 'a. 



4-2 



'c2.*d2+2.*a3.*b3 
'a2.*bl.*cl.*d2+2 
►c4.*d2-2.*al.*b4 
►a3.*b2.-c3.*d2-2 
'c2,*dl-2,*a4.*b3 
'a4.*b4.'cl.^'dl-2 
'c4 . 'dl-2 . 'a2 . 'b4 



'cl . 'd4-al . 'b2 . 'c4 
-a2 . *bl . 'c4 . ■d4+a3 . 'b4 . 'c4 . 'd4-a4 . 'b3 . ' 
b4 . ■'c2 . 'd4-a2 . 'b3 . ■c2 . 'd4+a3 . 'b2 . •c2 . 'd 
+a4 . *bl . 'c2 . *d4-a2 . 'b4 . 'cl . ■'d4+a4 . *b2 . ' 
b4 . *c4 . *d3-a4 . *b3 . *c4 . *d3+a2 . *b2 . *c3 . *d 
+a4 . *b4 . *c3 . *d3+al . *b4 . *c2 . *d3-a2 , *b3 . * 
b4 . *cl . *d3+a4 . *b2 . *cl . *d3; 

al-real (K(3iib2±nd(sM, 2'X-l, 2*X8 ( : , s) -1) ) ) ; 
a2-imag (K (3iib2±nd (sM, 2'X-l, 2*X8 ( : , s) -1) ) ) ; 
a3-real (K (sub2±nd (sM, 2'X, 2'X8 ( : , si -1) ) ) ; 
a4=imag (K (sub2ind (sM, 2'X, 2'X8 ( : , sj -1) ) ) ; 



.*a3.*b3.*cl.*dl-2. 
'c3.^'d3+a3.*b3.*c3. 

*d4-2. 
*d2-2. 
*d2-2 . 
*dl+2. 
■d4+2 . 
'd4+2 . 
d4+a2 . 'b2 . 'c3 . 
d4-al . 'b2 . 'c4 . 
d3+a3.*b2.*c2. 



*a4. *b4 
*c2.*d2 
*a3. *b4 
*c3.*d2 
*a4,*bl 
*c2,*dl 
*al. *b4 
•c3. 'dl 
'a3 . 'b3 
'c3 . 'd4 



bl.-*c3. 

1. *b3. ' 
b4.-*c2. 
b3.*cl. 
b2,*c2, 

03. *c4, 

02. *c3. 

33. *cl. 

4 . ■'b4 .' 
2.*bl.■ 
4.*bl.■ 



*dl-a3 . 'bl 
cl . 'd3 + 2 . ', 
'd4+2 . 'a4 . 
*d2+2.*al. 
-dl-2.*a2. 
*dl-2.'a3. 
*d4+2.*a2. 
*d4-2.*a3. 
c3.*d4+al. 
c4.*d3+a3. 
c2.*d3-a2. 



O ( (t-1) 'L"2 + k, 10, s) = ... 

-2 . *al . 'b3 . -cS. 'dl+2 . 'al . 'bl . *cl . 'dl-2 
bl . 'c3 . 'd4+a3 . -bS. 'c3 . 'd4 
•cl . *d4-a3 . 'bl . 'cl . 'd3-al . 'bl . 'cS. '■ 
c2 . 'd2+2 . 'a4 . •b4 . 'c2 . 'd2+2 . 'a2 . •b4 . '■ 



*cl . 'd4-, 
-al . 'b3 . 
.2.'b2. 'c: 
-2. *a4 . 
b2. 'c2. 
-2. *al. 
b3. *c4. 



*d2-2. *a2 



'dl-2 . 



■d4+2 . 



'd4+2 . 



b2 . 'cl . 
-a2 . 'bl 
b4 . 'c2 
+a4 .'bl 



b4.*cl 



d2+2.*a3.*bl. *c4 
b3.*c3.*d2-2,*a3 
■dl+2.*a3.*b4,*c2 
b2.*cl.*dl-2. *a4 
dl-2 . 'a4 . 'bl . 'c4 

d4-2 . 'a3 . 'b4 . 'c3 

b3 . 'c2 . ^d4-2 . 'a3 

d4-2 . 'a4 . 'bl . 'cl 

'C4.'d4+a3.'b4.'c4.'d4-a4.'b3.'c4 
'd4-a2 . 'bS . 'c2 . 'd4+a3 . •b2 . 'c2 . 'd4 

*c2 . *d4-a2 . 'b4 . 'cl . *d4+a4 . 'b2 . 'cl 
*d3-a4 . *b3 . *c4 . *d3+a2 . *b2 . *c3 . *d3 
. *c3 . *d3+al . *b4 . *c2 . *d3-a2 . *b3 . *c2 
*d3+a4.*b2.*cl.*d3; 



'b2 . *c4 . 'd4+2 . 'a4 
'd2+2 . 'a3 . *b3. 'c2 
'bl.'cl.'d2+2.*a3 



b3 . *cl . 'dl-2 

d3+a3.'b3. ■ c'i 
d2 ... 
*c4 . 'd4-2 

'cl. 'd2-2 

*c3.*d2-2 

*c4 . *dl+2 
'dl ... 
■b3 
■d4 ... 

. 'a4+2 
d4+a2 . 'b2 . 'c3 
d4-al.'b2.'c4 
d3+a3.*b2.*c2 



al-real (K (sub2ind {sM, 2*X-1 , 2*X9 ( : , s) -1) ) ) ; 
a2-iirag(K(sub2ind(sM,2*X-l,2*X9(:, s)-l) ) ) ; 
a3-real (K ( sub2ind (sM, 2'X, 2^'X9 ( : , s) -1) ) j ; 

a4=imag (K ( 3ub2ind (sM, 2'X, 2*X9 ( : , s) -1) ) ) ; 



bl = real (K (sub2ind (sM, 2'X8 ( 
b2 = imag (K (sub2iiid (sM, 2'X8 ( 
b3=real (K (sub2ind (sM, 2'X8 ( 
b4=imag (K (sub2ind (sM, 2'X8 ( 

cl-real (K (sub2ind (sM, 2*X9 ( 
c2-imag (K {3ub2ind (sM, 2*X9 ( 
c3-real (K (3ub2ind (sM, 2'X9 ( 
c4-imag (K (3ub2ind (sM, 2*X9 ( 

dl-real (K {3ub2ind (sM, 2*X10 ( 
d2^-imag (K<Biih2i nd (sM, 2'X1 ( 



-1, 2'X9 (: , 31-1) 1 ) ; 
-1, 2'X9 ( : , 3) -1) 1 ) ; 
;) ,2'X9(: , 3) -11 ) ) ; 
;) ,2'X9(: , 3l -11 ) 1 ; 

-1,2'XIO (:, si -11)1; 
0-1,2*X10 (;, s) -1) ) ) ; 
,2*X10{ :,s)-l) ) ) ; 
,2*X10{ :,s)-l) ) ) ; 

s)-l,2*X-l} ) ) ; 



bl = real (K { s^jb2ind (sM, 2'X9 ( : , 
b2=imag (K ( 3ub2ind (sM, 2'X9 ( : , 
b3=real (K ( sub2ind (sM, 2'X9 ( : , 
b4=imag (K ( siib2ind (sM, 2'X9 ( : , 

cl=real (K ( sub2ind IsM, 2*X10 I : 
c2-iinag (K (sub2ind {3M, 2*X10 { : 
c3-real (K (sub2ind {3M, 2*X10 { : 
c4-iiiiag (K [ sub2ind {3M, 2*X10 { : 

dl-real (K (sub2ind {3M, 2*X11 { : 
d2^i Tn^g (K ( snh2 -j nd < ■■gM, ? 'XI 1 { : 



= 1 -1, 2'XIO ( ; , s) -1) 1 1 ; 
3l -1, 2'XIO ( ; , s) -1) I 1 ; 
3l ,2'XIO ( : , 3l -II 1 ) ; 
3l ,2'XIO (: , 3| -II 1 ) ; 

.s)-l,2'Xll (: ,s) -1) 1 ) ; 
,s)-l,2*Xll (:,s)-l) ) ) ; 
,3) ,2'Xll (;, s) -1) ) ) ; 
,3) ,2'Xll (;, s) -1) ) ) ; 

,3)-l,2*X-l) ) ) ; 



Figure C.5: Programs used for the lattice calculation. 
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ol.*c3 

1. *b3. 
o4.*c2 
o3.*cl 

02. *c2 

03. *c4 
32.*c3 

03 . -cl 

4 . *b4 . 
2.'bl. 
4 . *bl . 



'd4+2. 



'd2+2. 



d3-real (K {3ub2ind (sM, 2*X11 ( : , s) ,2*X-1) ) ) ; 
d4-imag (K {3iib2ind (sM, 2*X11 ( : , s) ,2*X-1) ) ) ; 

0((t-l)*L'-2+k,ll,s)= ... 

-2 . *al . *b3 . *c3 . *dl+2 . *al . *bl . -cl . *dl-2 . *a3 . 
. *dl-a3 . -bl . *cl . *d4-al . *bl . *c3 . *d4+a3 . -b3 . *c3 . *d4 ... 

-al .*b3 . *cl . *d4-a3. *bl . *cl . -dS-al . *bl . *c3 . 
'cl,*d3+2 . *a2 . *b2 . *c2 . *d2+2 . 'a4 . *b4 . ■'c2 . *d2+2 . *a2 . *b4 . 

-2 . ■*a4 . 'b2 . *c4 . *d2 + 2 . *a2 . *b2 . ■'c4 . *d4 + 2 . *a4 
. ■*c2 . *d4-2 . *al . ■'bl . 'c2 . *d2+2 . *a3 . *h3 . • c2 
. *a.l . *b2 . *cl . ■d2-2 . '32 . *bl . *cl . •d2 + 2 . *a3 . 
. ■*c4 . 'd2+2 . 'a3 . *bl .'cl. 'd2-2 . 'al . *b4 .'c3. 
+2.-'a2.*b3.*c3.-*d2-2.-*a3.-b2.*c3.-*d2-2.-*a4. 
.*dl-2 .*a2 . *bl . *c2 . 'dl+2 . 'a3 . *b4 .'cl. 'dl-2 . 'a4 . *b3 .'cl. 

-2 . *a2 . *b2 . *cl . *dl-2 . *a4 . *b4 . *cl . *dl-2 . *al . 
. *dl-2 . *a3 . *b2 . *c4 . *dl-2 . *a4 . *bl . *c4 . *dl-2 . *a2 . *b4 . *c3 . 

+2 . *a4 . *b2 . -c3 . *dl-2 . *al . *bl . -c4 . *d4+2 . *a3 . 
*d4+2 . *a2 . *bl . *c3 . *d4-2 , *a3 . *b4 . *c3 . *d4+2 , *a4 . *b3 . *c3 . 

-2 . *al . *b3 . *c2 . *d4-2 . *a3 . *bl . *c2 . *d4-2 . *al . 
*d4-2 . *a3 . *b2 . "'cl. *d4-2 . *a4 . *bl . *cl. *d4-al . *b2 . *c4 , *d4 

-a2 . 'bl . *c4 . *d4+a3. *b4 . *c4 . *d4-a4 , *b3 . *c4 
■c3. *d4+al . *b4 .'cl. *d4-a2 . *b3 . *c2 . *d4+a3 . *b2 . *c2 , *d4 

+a4 . *bl . *c2 . *d4-a2 . *b4 .'cl. ■*d4+a4 , *b2 . *cl 
■c4 . *d3 + a3 . *b4 . •c4 . 'd3-a4 . *b3 . ■'c4 . 'd3+a2 . *b2 . *c3. *d3 

+ a4 . *b4 . 'c3 . 'dS + al . 'b4 . *c2 . ■'d3-a2 . *b3 . *c2 
'c2 . *d3-a2 . *b4 . 'cl . 'd3+a4 . *b2 . 'cl . 'd3; 

al-real (K (sub2ind (sM, 2'X-l, 2'Xll ( : , sj -1 j ) ) ; 
a2-imag (K [sub2±nd(sM, 2*X-1, 2'Xll ( : , s) -1) ) ) ; 
a3-real (K (sub2ind (sM, 2*X, 2*X11 [ :,s)-l) ) ) ; 
a4-imag (K (sub2ind (sM, 2*X, 2*X11 ( : , s) -1) ) ) ; 



b3.-cl.*dl-2. 

d3+a3.*b3.*c3. 

c4. -d2 
*b4 . ■'C4 
*d2 

'd2-2. 
•d2-2. 
*dl+2 . 
*d4+2. 
*d4+2. 
d4+a2,*b2.*c3. 
d4-al.*b2.*c4. 
d3+a3.*b2.*c2. 



% performs calculations on saved configs 
:;lear all 

global Psi D L N Xc Xs Xp Xm Xop.=; 

inf ile=' /home/wytham/r johnson/myf ree/ ops7/ops7b611 6' , 
Dutf ile=' /home/wytham/r johnson/ lattice/Htest3ops7p' ; 



bl-real (K (sub2ind (sM, 2*X11 ( ; , 
b2-imag (K {siib2ind (sM, 2*X11 ( ; , 
b3-real (K {siib2ind (sM, 2*X11 ( ; , 
b4-imag (K {sub2ind (sM, 2*X11 ( ; , 

cl=real (K (sub2ind (sM, 2*X12 ( : 
c2 = imag {K (subzi.Td (sM, 2*X12 ( : 
c3=real (K (sub2ind (sM, 2*X12 ( : 
c4-iraag (K [sub2ind(sM, 2*X12 ( : 



s)-l,2*X12(:,s)-l))t; 
st-l,2*X12 { :,s)-l) ) t ; 
St ,2*X12 (:,s)-l) ) ) ; 
s) ,2*X12 (: ,s)-l) ) ) ; 

sj -1, 2*X7 (: , sj -Ij ) ) ; 
sj -1, 2-'X7 ( : , sj -Ij ) 1 ; 
sj ,2'X7 ( : , s)-l) 1 j ; 
sj ,2-'X7 [ :, s)-l) ) ) ; 



dl-real (K {sub2ind (sM, 2*X7 ( : , s) -1,2*X-1) ) ) ; 
d2-imag (K (sub2ind (sM, 2*X7 ( : , s) -1,2*X-1) } ) ; 
d3-real (K (sub2ind (sM, 2*X7 ( : , s) , 2*X-1J ) ) ; 
d4 = imag (K (sub2ind (sM, 2*X7 ( : , s) , 2*X-1J ) ) ; 



( (t-1 j *L"2+k, 12, sj = ... 

-2 . 'al . *b3 . *c3 . ■'dl + 2 . "'al . *bl . *cl . 
[bl . *c3. *dl-a3 . *bl . *cl . *d4-al . *bl . 'cS. *d4+a3 . *b3 . 'cS 

-al . *b3 . *cl . ■'d4-a3. 'bl . *cl . 'dS-al . *bl . *c3 . 
b3 . *cl . *d3+2 . *a2 . *b2 . ■*c2 . *d2+2 . *a4 . *b4 . *c2 . *d2+2 . *a2 . *b4 . 

-2 . ■*a4 . *b2 . *c4 . *d2+2 . *a2 . *b2 . *c4 . *d4+2 . *a4 
*c2 . *d4+2 . *a4 . *b2 . *c2 . *d4-2 . *al . *bl . *c2 . *d2+2 . *a3 . *b3 . *c2 
-2 . *al . *b2 . *cl . *d2-2 . *a2 . *bl . *cl . *d2+2 . *a3 
*cl . *d2+2 . *al . *b3 . *c4 . *d2+2 . *a3 . *bl . *c4 . *d2-2 . *al . *b4 . *c3 
+2 . *a2 . *b3 . *-c3 . *d2-2 . *a3 . *b2 . *c3 . *d2-2 . *a4 
*c?.*d1-?.*a2.*b1 .*c?.*d1+?.*a3.*b4.*c?.*dl-?.*a4.*b3.*c2 



*b2.*cl.*dl- 
'(II 



dl-2.*a3.*b3.*cl.*dl-2.*a3. 



d3+a3.*b3.*c3.*d3-i 
c4.*d2 ... 
*b4.*c4.*d4-2.*a2. 
*d2 ... 

*b4.*cl.*d2-2.*a4. 
*d2 ... 

*bl.*c3.*d2-2.*al. 
*d1 ... 

*b4.*c4.-'dl + 2.*a2. 



aval ( [ ' delete ' 
% eval ( ['load ' 



outfile num2str(P) 
infile ' 1 .mat' ] ) 



Hsweeps-2000; 

alph-[.l .25 .26 .15 .2]; 

la=le_Tqtri (alphj ; 

% set up operators 

B-1; 

L-16; 

D-3; 

H-2; 

ops- [2 3 4 5 6]; % 
lops=length(ops) ; 



;<p=repmat (Xc, [1 D) ) ; 

:<p [ l)-sub2ind( [L L L] , mod (Xs ( : , 1 j , Lj +1 , Xs [ : , 2) , Xs ( : , 3 j ) ; 

Kp(:,2)-sub2ind( [L L L] , Xs ( : , 1 ) , mod (Xs ( : , 2) , L) +1, Xs ( : , 3) ) ; 

Kp{:,3)-sub2ind( [L L L] , Xs ( : , 1) , Xs { : , 2) ,mod (Xs { : , 3) , L) +1) ; 



repmat (Xc, [1 D] ) ; 

: , 1) =sub2ind ( [L 1. L] , 

: , 2) =sub2ind ( [L 1^ L] , 

3)=sub2ind( [L L L] , 



; , 11 -2, Lj +l,Xs ( : , 21 ,Xs ( : , 31 j ; 
,mod{Xs ( ; , 2} -2, L) +l,Xs ( : , 3| j ; 
, Xs ( ; , 2) ,mod (Xs ( : , 3) -2, Lj +1 j ; 



3pp0=ze: 
3pp2=ze: 
3pp4=ze: 



bnO=7.eros (NsweeDs/ski p, 1 ods. 1 a I ; 

irppO^zeros iNsweeps/skip. lops . j.ops. L/2 + 1, la) ; 



OS {Nsweeps/skip, la) ; 

s(L,lops,lat; 
s (L, lops, la) ; 
s (L, lops, la); 



bl.-c4. *d3 + a3.*b4.*c4. -d3-a4 .*b3.-*c4. -d3+a2.*b2.*c3.*d3 . . . 

+a4 . *b4 . 'c3 . ■'d3+al . 'b4 . *c2 . ■'d3-a2 . *b3 . *c2 . *d3+a3. *b2 . *c2 . *d3+a 
4 . *bl . *c2 . *d3-a2 . *b4 . ■'cl . *d3+a4 . *b2 . *cl . *d3; 

end % for s 

kplaq(h,at-(mean(reaHdiag(K) ) J -l-4*alph (at "2-28*alph (a) "4) /8/alphla 



end % for t 
for t-l:12 

nc(t, :,a)-sqrt (mean(0{ :,1, : t ."2) . /mean (0 ( : , t , :) .^2) J ; 

% take out the normalization of the loops so it 

% for t-l:12 
% for 1-1: lops 

% 0(:,t,l)-0(:,t,l) .*nc(t,l,a) 

% end 

% end 



! done later . . 



% 0-re shape (O, L-2, L, 12, lops) ; 

% pp=squeeze (mean (Oj ) ; 

pp-pp+0; 

if pk=-l 

ppp=pp; 

pp=sqjcezc (^-ican (reshape {pp, L''2, L, 12 , lops) ) ) ; 
eval(['save ' outfile 'pla' niam2str(a) 'h' num2str(h) '.mat alph i 
■ flag Xops Xc' ] ) 
pp-ppp; 

save calcl6ops2.dat h -ascii 



p±iippi ztJl'_l^: (U, lop'j, J.ti| ; 
oinC2=zeros (B, lops, lops, L/2+1 , la) ; 
jackC2=zeros (B, lops, L/2+1, laj ; 
jackE2-zeros (B, lops, L/2+1, laj ; 
jackM2-zeros (B, lops, floor (L/3j , la) ; 

pp4=zeros (Nsweeps/ skip, lops, la) ; 

orpp4=zeros (Nsweeps/skip, lops, lops, L/2 + 1, la) ; 
binpp4=zeros (B, lops, la) ; 
binC4-zeros IB, lops, lops, L/2+1, la) ; 
jackC4-zeros (B, lops, L/2+1, la) ; 
jackE4-zeros (B, lops, L/2+1, la) ; 
jackM4-zeros (B, lops, floor (L/3) , la) ; 

%0++/4++ operator 

orpp04=zeros (Ns weeps/ skip, lops, lops, L/2+1, la) 
|binC04=zeros (B, lops, lops, L/2+1, la) ; 

%0++/2++ operator 

orpp02=zeros (Ns weeps/ skip, lops, lops, L/2+1, la) 
|binC02-zeros (B, lops, lops, L/2+1, la) ; 

%2++/4++ operator 
■orpp24-zeros (Wswocps 'skip, lops, lops, L/2 + 1, la) 
|binC24=zeros (B, lops, lops, L/2+1, la) ; 



for h=l : Ns weeps/ skip 
disp (h) 

^for a-l:la 



%4++ operate 



a' num2str(a) 
num2str (a) ' h 



h' num2str { (h-1) *skip+l } 
num2str ( (h-1) *skip+l| ' 



eval (['save ' outfile ' 
flag Xops Xc' ] j 

% save /home/wytham/r johnson/myf: 
pp-zeros (length (Xc) , 12, lops) ; 



num2str{a) ' h' num2str{h) ' .mat alph ops pp j 
/new4 recover .mat 



3nd % for h 

avgplac[=mean (kplaq) ; 

tocc=etime (clock, ticcj ; 

% clear unnecesary variables 

%clear Psi ppO corppO ppl corppl pp2 corpp2 pp3 corpp3 pp4 corpp4 
|clear Psi 

Ll[['save ' outfile 'tot. mat ticc tocc avgplaq' ] ) 



pp^squeeze (mean ( reshape (pp, L''2,L, 12, lops) ) ) ; 
%for t-l:12 

% nc(t, :,a)-sqrt (mean (pp ( : , 1, : ) . "2) . /mean (pp (mod ( [1:L],L)+1,- 

for t=l:12 

for 1=1 : lops 

pp(:,t,l)-pp(:,t,l) .*nc{t,l,a) ; 



r-ops (1) ; 

ra^rnnnd (r*Hqrt ( 3) ) ; 



Figure C.6: Programs used for the lattice calculation. 
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2-atan {r/ra) ] , 



%th-[atan{r/ra) pi-atan (r/ra) 3*pi/2 -atan{r/ra) pi/2 pi+atan (r/ 
% 3*pi/2+atan (r/ra) pi/2-atan (r/ra) pi pi/2+atan (r/ra) 3*pi/ 
th-pi.*[l/6 5/6 3/2 11/6 1/2 7/6 10/6 1/3 1 2/3 4/3] 

J-0; ph-exp(i*J.*th); 
OppO (t, l,a) - ( ... %niean { ... 

pp(t,l,l)+pp{t,2,l)+pp(t,3,l)+pp(t,7,l)+pp{t 



)/12; 
,1) +0(: 



■,1> )/3; 
.,1) )/6; 



■,1} )/3; 
)/6; 



% pp(t, 
% 

, 5,l)+0(:, t, 6,1) ) /6; 



l)+pp(t,9,l) )/ 

l)+pp(t,S,l)+pp(t,6,l)+pp{t,10,l)+pp(t,ll,l)+PP(t,12,l 
t,l,l)+0(:,t,2,lj+0{:,t,3,l) )/3; %+0{: 



■-If )/3; 
)/6; 



J-1; ph-exp(i-J.-th); 
Oppl(t,l,a)-mean( ... 

ph(l)*0(:,t,l,l)+ph(2)*0(:,t,2,lj+ph(3t*0(:, 

ph{4)*0(;,t,4,l)+ph(5)*0(:,t,5,l)+ph(6)*0{:. 



J=2; ph=exp (i*J. *th) ; 
Opp2 (t,l, a]= ( ... %mean ( 

ph(l) *pp(t, 1,1) +ph (2) -pp (t, 2, 1) +ph(3j -pp(t, 3, 1) + ... 

ph(7) *pp(t, 7,1) +ph (8) -pp (t, 8, 1) +ph(9j -pp(t, 9, 1) j /6; 
% ph{4) *pp[t, 4, 1) +ph(5j -pp(t,5, l)+ph(6) *pp(t, 6,1) + ... 
% ph{10) *pp(t, 10,1) +ph (11) *pp(t, 11,1) +ph (12) *pp(t, 12,1) ) /12; 
% ph[l) ■'0(:, t, l,l)+ph(2) *0(: ,t,2,l) +ph(3) *0l 

% ph{4)*0(;,t,4,l)+ph(5)*0(:,t,5,l)+ph(6)*0{:, 

% J-3; ph-exp (i* J. *th) ; 

% Opp3(t,l,a)=mean( ... 

% ph(l) ■'0(:, t, 1, l)+ph(2) *0(: , t, 2, 1) +ph(3) '0 [ : , " 

% pli{4) *0(:, t,4,l)+ph(5) *0 ( : , t , 5, 1) +ph (6) *0{ :, 



J-4; ph-exp(i*J.*th); 

Opp4 (t, l,a) - ( ... %mean{ . . . 

ph(l) -pp(t, l,l)+ph{2) *pp(t, 2, 1) +ph(3) *pp(t,3, 1) + ... 

ph(7) *pp(t, 7,l)+ph (S) ---pp (t, 8, 1) +ph(9) *pp(t, 9, 1) ) /€; 
% ph{4) *pp(t, 4, 1) +ph(5) *pp(t, 5, 1) +ph(6) -'pp(t, 6, 1) + ... 
% ph{10)*pp(t,10,l) +ph (11) *pp(t, 11, 1) +ph (12) *pp(t, 12, 1) ) /12; 
% ph[l)-0(:,t,l, l)+ph(2) *0(: , t, 2, 1) +ph(3) -O ( : , ■ 

% ph{4)*0(;,t,4,l)+ph(5)*0(:,t,5,l)+ph(6)*0{:. 



C^zeros (lops, lops) ; 
for t-0:L/2 

for 1=1: lops 

for k=l:lops 

C(l,k) -mean (OppO ( : , 1, a) . *OppO (mod{ [l:L]+t-l, L)+l,k,a) ) ; 



%ppO=reshape (ppO, [Nsweeps/B B lops la]); 

rppO=reshape (corppO , [tJsweeps/B B lops lops L/2+1 la]); 



%pp2=reshape (pp2, [Nsweeps/B B lops la]}; 

%corpp2=reshape (corpp2, [Nsweeps/B B lops lops L/2+1 la]); 
notions 

binC2=squeeze (mean (corpp2 ) ) ; 
binpp2=squeeze (mean (pp2 ) ) ; 

%pp4=reshape (pp4, [Nsweeps/B B lops la]); 

%Corpp4=reshape (corpp4, [Nsweeps/B B lops lops L/2+1 la]); 
notions 

binC4=3queeze (mean (corpp4 ) ) ; 
binpp4=scpaeeze (mean (pp4 ) ) ; 



binC02=scpaeeze [mean (corpp02) ) ; 
binC04=squeeze (mean (corpp04) ) ; 
binC24=squeeze {mean (corpp24) ) ; 



% meanplaqO (k, a) =squc-c' 

% stdplaqO (k, a) =squeezf 

% meanplaq2 (k, a) =squee: 

% stdplaq2 (k, a) ^squeeze (std [binpp2 (k, 

% meanplaq4 (k, a) ^squeeze (mean (binpp4 (k, a) 

% stdplaq4 (k, a) =squeeze (std (t>inpp4 (k, a) ) ) 

meanplaqC (k, a) ^bi-ppC (k 
stdplaq(j (k , .5) t a (opC { 



binppO (k, a) 
.nppO (k,a) ) ) 
n (binpp2 [k,a) 



meanplaqO (k, a ) -oi:".pp2 { 
stdplaq2 (k, a) =std (pp2 ( 
meanplaq4 (k, a) =binpp4 (k 
stdplaq4 {k,a) -std{pp4 ( 



) ) ; 



: subtracting here,. 



prevacO = binCO; 
prevac4 — binC4; 

for b-l:B 

for 1=1: lops 

for k^'l : lops 
binCO (1, k, 

for b-l:B 

for 1=1: lops 

for k=l : lops 
binC4 (l,k. 



, a) =binCO (1, k, : , a) - binppO (1, a) . 'binppO (k, a) ; 



a)-binC4 (l,k, :,a) - conj (binpp4 (l,a) ) . *binpp4 (k, a) ; 



C=zeros ( lops, lops) ; 
for t-0:L/2 

for 1=1: lops 

for k=l:lops 
C(l,k)-mea 



(real (Opp2 { :,l,a) .*conj (Opp2 {mod ( [1 : L] +t-l, L)+l,k, 



pp2{h, :,a)-mean{0pp2(:, : , a) ) ; 
corpp2 {h, :, : , t+1, a) = (C+C . ' ) /2; 



C^zeros ( lops , lops ) ; 
for t-0:L/2 

for 1=1 : lops 

for k^l : lops 
C(l,k)-mea 



for b-l:B 

for 1=1 : lops 

for k=l : lops 
binC02 (l,k, 
binC02 (l,k. 



a) ] 



(real (Opp4 { :,l,a) .'conj (Opp4 {mod ( [1 : L] +t-l, L)+l,k, 



(OppO ( :,l,a) .*conj (Opp4 (mod ( [ 1 : L] +t-l, L) +l,k,a) ) ) ; 



pp4 {h, ;,a)-mean{0pp4 (:, : , a) ) ; 
corpp4 {h, :, : , t+1, a) = (C+C . ' ) /2; 

C^zeros ( lops, lops ) ; 
for t-0:L/2 

for l=l:lops 
for k=l:lopi 
C(l,k)-m. 

corpp04 (h, :, : , t+1, a) - {C+C . ' ) /2; 

C^zeros (lops, lops) ; 
for t-0:L/2 

for 1=1 : lops 

for k=l : lops 

C (l,k)=mean(OppO ( : ,l,a) . *con j (Opp2 (mod ( [ 1 : L] +t-l , L) +l,k,a) ) ) , 

corpp02 (h, :, : , t +1, a) - [C+C . ' ) /2; 

C^zeros (lops, lops) ; 
for t-0:L/2 

for l=l:lops 
for k=l:lop: 
C(l,k)-mi 

corpp24 (h, :, : , t+1, a) - {C+C . ' ) /2; 

% save -/myfree/tmpsave .ma. 

id 

-append 



for b=l:B 

for 1=1 : lops 

for k=l : lops 
binC24 (l,k, 
binC24 (l,k, 

l,a)); ^ 



■ b=l: 
x=l ;B 



a)-binC02(l, 
a)-binC02 (1, 



,a)-binC24(l, 
,a)-binC24(l, 



, a) - conj (binpp4 (k, a) ) . 'binppO (1 , a } ; 
,a) ./sqrt (binC4 (l,k, l,a) ) . /sqrt {bi.iCO (1, 



,a) - con j (binpp2 (k, i 
,a) ./sqrt (binC2(l,k, 



a) - conj {binpp4 (k, a) ) . *binpp2 (1, a t ; 
a) ./sqrt (binC4 (l,k, l,a) ) ./sqrt{binC2 (1, 



jackCO (b, 
jackCO (b, 
jackEO (b, 
jackMO (b, 

jackC2 (b, 
jackC2 (b, 
jackE2 (b, 
jackM2 (b, 

jackC4 (b, 
jackC4 (b, 
jackE4 (b, 
jackM4 (b, 



a) ^squeeze (mean (binCO (x, k, k, : , a) ) ) ; 
a)-jackCO{b,k, :,a) ./ jackCO (b, k, 1, a) ; 
a) -squeeze (std(binCO (x,k,k, : , a) ) ) ; 
a)-log{jackC0(b,k,l:L/3,a) . /jackCO (b, k, 1+ [1 : L/ 3 

a) -squeeze (mean (binC2 (x, k, k, : , a) ) ) ; 

a) = jackC2 (b, k, ; , a) . / jackC2 (b, k, 1, a) ; 

a) =squeeze (std(binC2(x,k,k, :,a))); 

a) =log ( iackC2 (b,k,l;L/3,a)./ jackC2 (b,k,l+[l;L/3 

a) =squeeze (mean (binC4 (x, k, k, : , a) ) ) ; 
a) = jackC4 (b, k, ; , a) . / jackC4 (b, k, 1, a) ; 
a) -squeeze (std(binC4 (x, k, k, : , a) ) ) ; 
a)-log{jackC4(b, k,l:L/3,a) . / jackC4 (b, k, 1+ [1 : L/ 3 



jackMO (b, k, a) -fmins { ' f itmass' , 1 ,[],[] , - - - 

squeeze [ jackCO (b, k, : , a) ) , squeeze ( jackEO (b, k, : , a) ) ) ; 



(Opp2 ( : ,l,a) .-con j (Opp4 (mod ( [ 1 : L] +t-l , L) +l,k,a) ) ) ; 



! calcl6err.dat h 



% save 
3nd % for h 

this is where herr 



% 








% 


meanCO (k 


: , a) -mean (squeeze ( jackCO ( 


,k. :,a) ) ) ; 


% 


stdCO (k. 


, a) -sqrt (B-1 ) *3td ( squeeze 


( ;ackCO ( : , k 




mOpp (k, : 


a) -mean ( jackMO ( : , k, : , a) ) ; 






sCpp (k, : 


a) =sqrt (B-1) *std ( jackMO ( ; 


k, : ,a) , 1) ; 




moanC2 (k 


: , a) =mean (squeeze ( jackC2 ( 


,k, : , a) ) ) ; 




stdC2 (k, : 


a) =sqrt (B-l)*std(squeeze(jackC2 (:,k. 




m2pp (k, : 


a) =mean ( jackM2 ( ; , k, : , a) ) ; 






s2pp (k, : 


a) =sqrt (B-1 ) 'std ( jackM2 ( ; 


k, : ,a) , 1) ; 




meanC4 (k 


: , a) =mean (squeeze ( jackC4 ( 


,k, a) ) ) ; 


% 


stdC4 (k. 


, a) -sqrt (B-1 ) *std (squeeze 


( jackC4 (: ,k 


% 


m4pp {k, ; 


a) -mean { jackM4 ( : , k, : , a) ) ; 




% 


s4pp{k, : 


a) -sqrt {B-1) -std { jackM4 ( : 


k, ;,a),l); 



/home/wytham/ r johnson/myfree/newSrecover .mat 

he used to reninve had cnnfigs, T think 



meanCO (k, : ,a) -squeeze (binCO (k, k, ; , a) ) ; 
meanCO (k, : ,a) -meanCO {k, ; , a) . /meanCO {k, 1, a) ; 
mOpp(k, : ,a)-1ng(TnearirO (k, 1 :T,/3,a) . /meanCO { k , 1 + 



Figure C.7: Programs used for the lattice calculation. 
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meanC2 (k, : , a) -squeeze (binC2 (k, k, j , a) ) ; 

nieanC2 (k, : ,a) -meanG2 (k, ; , a) . /meanC2 (k, 1, a) ; 

m2pp(k, : ,a)-log(meanC2 {k,l;L/3,a) . /meanC2 (k,l+[l:L/3] , 

meanC4 (k, : , a) -squeeze (binC4 (k, k, : , a) ) ; 

meanC4 (k, : , a) -meanC4 [k, : , a) . /meanC4 (k, 1, a) ; 

m4pp (k, : , a) -log (meanC4 {k, 1 : L/3, a) . /meanC4 (k, 1+ [1 : L/3 ] , 



LnCD4 (k, : , a) =3queeze (blnC04 (k, k, 

LnC02 (k, ; , a) =3queeze (binCD2 (k, k, 

meanC02 (k, :,a)-meanC02 (k, 
LnC24 (k, : , a) -squeeze (binC24 (k,k, 

meanC24 (k, : , a) -meanC24 {k, : , a) - /meanC24 (k, 1, 



:nea_TC04 (k, 1, a) ; 
/meanC02 (k, 1, a) ; 



end 



end % a 

avgplaq-mean (kplaq) 
%clear Psi ppO co 



rppO ppl corppl pp2 corpp2 pp3 corpp3 pp4 corpp4 



outfile num2str(P) 



function V-update J (sites, mu, R) ; 
% updates V<-U given R 

global Psi Parameters Xc Xs Xp Xm 

H-Parameters .N; 

■ml (Psi ( (mu-1) *N^2+ [l:N-^2] , sites) , 

k-sqrt (sum(Z.'2) ) ; 
u-Z./repmat(k, [N-'2 1]); 



alpha— 2*Paramete: 



i.Beta 



*k./N; 



= zeros (1, length (sites) ) ; 
■ 1=1 : length (sites) 

% generate aO s.t. P (aO) - exp{alpha aO) (l-a0'2)'.5, {-1<^ 

r4=l; dl=2; 

while (r4-2 > (l-dl/21) 

rl-rand; r2-rand; r3-rand; r4-rand; 

xl — log(rl) /alpha (1) ; x2 — log {r2 ) /alpha (1) ; 

c-(cos (2*pi*r3) ) "2; 

a-xl*c; % b-xl-a; for another ind aO 

dl-x2+a; % d2-x2+b; 



end 

rl-2*rand(size (aO) )-l; r2-rand(s 

r3-3qrt (1-aO . '2) ; 

a3=rl .*r3; 

r4=3qrt (1-rl . "2) ; 

al = r3 . *-r4 . V-o-, C^'x i . ' -Z) ; 



function V-ovrlaxJ (sites, mu, R) ; 

% performs over relax update V<-U for SU{2) 

global Psi Parameters Xc Xs Xp Xm 

H-Parameters .N; 

■sqrt (sum(R.^2) ) ; 
|u-R./repmat(k, [N~2 1]); 

i3 (u,m4 (Psi( [mu-l)*N'-2+[l:N"2] .sites) ,u) ) ; 



Figure C.8: Programs used for the lattice calculation. 
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V = 2 jt p a + ij—^) (1 -exp(-f p)) 
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